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The report documents techniques and a FORTRAN 77 computer 
program for the identification of flutter-mode frequencies and 
damping ratios in near-real-time. Recent advances in recursive 
identification algorithms and analysis are applied to the flutter 
test problem. The accuracy and convergence of the algorithms are 
predicted analytically and substantiated with results from simu- 
lated and flight data. The results show promise for monitoring 
aircraft flutter characteristics in real-time with a high degree 
of reliability. 

The work has been performed under NASA Contract NAS4-2955. 
Technical direction and discussions with NASA technical monitor 
Mr. Glenn Gilyard and with Mr. Richard Maine are gratefully ac- 
knowledged. / 
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SECTION 1 


INTRODUCTION 


Real-time monitoring of flutter parameters can significantly 
improve safety of flight tests for new and modified aircraft, 
reduce flight test time, and aid envelope expansion. The moni- 
toring system must be reliable and accurate and should require 
minimum inputs from the operator. To aid flutter clearance of 
many aircraft, it is desirable to monitor flutter characteristics 
with randomly forced (turbulence excitation) as well as known 
forcing inputs. The most success is likely to be achieved with 
the use of advanced system identification techniques. For wide- 
spread use, the computation time should be suitable for a real- 
time implementation and a standard computer language like the 
ANSI FORTRAN-77 is desirable. 

The goal of this effort is to develop a computer program-- 
based on an identification algorithm — with excellent accuracy and 
convergence performance in estimating flutter mode characteristics 
of test aircraft from real data. 

1.1 SUMMARY OF APPROACH 

A summary of the approach used to develop and document the 
near-real-time flutter analysis program is shown in Figure 1-1. 

The performance of two different algorithms was evaluated 
with a high-level command-driven identification and control pro- 
gram called MATRIX^ [1]. This capability enabled a thorough 
yet efficient evaluation of the algorithms before implementation 
and testing of the Near-Real-Time Flutter Analysis (NRTFA) pro- 
gram in FORTRAN-77 source code. 
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SOFTWARE SPECIFICATION & INTEGRATION 1 

• 
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• 
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DOCUMENTATION 


Figure 1-1 Real-Time Flutter Analysis Approach 
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1.2 RESULTS SUMMARY 


The selected techniques are based on recursive prediction 
error methods (RPEM) , which can emulate a maximum likelihood, 
instrumental variables and least squares. 

The Cramer-Rao estimate error accuracy bounds for a two- 
active mode model predict that the critical parameter, damping, 
can be estimated quite satisfactorily after approximately 500 
data points. Actual simulations indicate that the convergence 
transient from no a priori information lasts slightly longer, but 
that for long data records of 2500 data points, damping ratios 
can be estimated very accurately (<.05% error). Modal parameters 
can be estimated from unknown inputs (turbulence excitation) with 
sufficient accuracy when the signal-to-noise ratio is approximately 
ten or greater. For lower turbulence levels almost any known input, 
even at very low levels, significantly improves the estimation 
accuracy . 

The delivered algorithm tracks the flutter mode of the third 
DAST flight quite successfully, indicating where the damping 
trend breaks toward zero at the onset of a pulse response decay, 
prior to the flutter condition. Quantitative comparisons of the 
convergence, robustness and accuracy of the recursive identifica- 
tion algorithm analysis and simulation results show the near— real- 
time flutter analysis (NRTFA) program to be highly reliable. 

1.3 REPORT ORGANIZATION 

Candidate algorithms and the associated model forms most 
suitable for use with them in flutter testing are described in 
Section 2. The reliability, robustness and accuracy of these 
flutter mode estimates are evaluated in Section 3, supported by 
results on simulated and flight data in Sections 4 and 5. The 
program architecture, configured with an interactive driver which 
calls the appropriate identification subroutines, with an instal- 
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lation guide is documented in Section 6. The memory and CPU 
requirements for the flutter analysis program as it is now con- 
figured are also given. Appendix A serves as a brief User's 
Guide giving an example interactive setup session, with guide- 
lines for parameter setup options, and a description of the modal 
parameter estimate output. Appendix B gives the FORTRAN listings 
of the principal routines of the NKTFA program. 
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SECTION 2 


REQUIREMENTS FOR REAL-TIME FLUTTER ANALYSIS 


The development of a system identification software 
package to monitor aircraft flutter characteristics is 
influenced by the following considerations: 

1. The basic requirements in real-time flutter analysis 
involves estimation of damping ratios and changes 

in natural frequencies of aeroelastic modes. Mode 
shapes are not required but may be determined if 
they improve damping ratio and natural frequency 
estimates. Physical parameters are usually not 
estimated in real-time. 

2. Several flutter modes must be tracked simultaneously. 
The modes can be closely spaced. For successful 
real-time implementation, computational efficiency 

is very important. 

3. The real-time software must be robust and reliable 
because its failure could cause delay or termina- 
tion of a flight test. 

4. Since the flight test engineer is generally busy 
with several tasks during the flight test, the 
software should require minimum commands during the 
test. Set-up parameters, if any, should be entered 
prior to the flight test, and the algorithm should 
be able to restart without entering any startup 
values . 

5. Because of the damping ratios of flutter modes are 
required with high precision and the data is often 
noisy, advanced identification techniques are 
highly desirable whose convergence properties and 
accuracy performance can be analyzed. 

6. The estimated parameters must be updated at regular 
intervals to continuously track any changes in 
flutter characteristics. Thus, a recursive tech- 
nique is desirable. 
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The following characteristics of a real-time flutter- 
monitoring software package determine the extent to which 
the above requirements are met. 

1. model representations, 

2. identification procedures and robustness 
characteristics, and 

3. numerical techniques, program architecture and 
input/output structure. 

These issues are interrelated since identification 
procedures and numerical techniques depend on the model form. 
The model form affects the shape of the criterion function 
we seek to minimize with respect to the parameters, as well 
as the type of equations which must be propagated to find 
the gradients. This section explains why certain model 
forms and algorithms are preferred for real-time flutter 
analysis. The identification algorithms are broadly grouped 
as recursive and block recursive algorithms. Robustness 
characteristics, numerical techniques and a suitable program 
architecture are discussed in subsequent sections. 

2 . 1 Model Forms 

Three model forms can be used: 

1. transfer function or continuous autoregressive 
moving average (ARMA), 

2. sample-data transfer function or discrete ARMA, 

3. state variables in continuous or discrete forms. 

2.1.1 Transfer Function 

The q x 1 input u to p x 1 output y transfer 
function is 

= T(s) , (2.1) 

u( s ) v v ' 
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which may be written as 


T(s) 


N(s) 

D(s) 


E 


i=l 


R. 

i 


(s 4 X.) 


( 2 . 2 ) 

(2.3) 


N(s) is a matrix polynomial and D(s) is a scalar poly- 
nomial in s . R. are pxq matrices and X. are the 

i l 

pole functions. 

For scalar input and output, the transfer function is 
also written as 

m 

T(s) = K £ (s 4 z.) ^ (s 4 x.) , (2.4) 

i=l i=l 

where z^ are the zeros and K is the gain. A continuous 
autoregressive moving- average form is 
n m 

y (n) a y^- 1 ' =£ B. u'"- 1 ' , (2.5) 

i=l 1 i=o 

where y ^ is the jth derivative of y . 


2.1.2 Sampled Data Transfer Function 

This transfer function and ARMA representation is 
similar to the continuous transfer function except a z-trans- 
form is used in place of the s-transform. Thus 

y(z -1 ) = T( z -1 ) u(z -1 ) (2.6) 


where 

z '\ = y k -i (2 - 7 

This model can be written in any of the forms corresponding 
to Equations (2. 3) to (2.5). 
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The representation maps the imaginary axis into the 
unit circle. Lightly damped oscillatory continuous modes 
will be close to the unit circle in the ARMA formulation, 
depending on how far the continuous roots are into the left 
half plane and the sampling time. The continuous damped 
frequency and sample time determine the angular orientation 
of the poles in the z-plane , 

0 = Tan -1 cj (1 - E, 2 )^ t . 

The different discrete polynomial models for both known 
and unknown inputs and their state-space realizations are 
discussed more fully in the next subsection on recursive 
algorithms, where they are particularly attractive. 

2.1.3 State Variable Models 

Defining the state vector as x , a state variable 
model is written as 

x = Fx + Gu + Tw , (2.8) 

y = Hx + v 

w is the process noise vector and v is the measurement 
noise vector. F, G and H matrices can take various forms 
depending on the particular selection of the state vector. 

In many offline applications, the state vector is selected 
to represent a set of physical quantities (e.g., body rates, 
deflection at a certain point on the wing). Such represen- 
tations give rise to physical parameters, such as stability 
and control coefficients. 

In real-time flutter analysis, the natural frequencies 
and damping ratios of aeroelastic modes are of primary 
interest. Since a state variable model [F, G, H] has the 
same input-output behavior as a [TFT - , TG, HT - ^] model if 
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T is a nonsingular matrix, the following representation 
appears useful (for single-input systems): 



( 2 . 10 ) 


0 

1 



o 

0) 

2 

-co 

n . 
i 

- 25 . y 

^i H n . 
l 

> 

or 

-CO 

a 


and H is a general matrix with potentially all nonzero 
elements. With multi-input systems, the remaining part of 
matrix G is general and unknown. 1 is a general matrix. 


The turbulence excited structural vibrations are of 
special interest. In this case G = 0 . For the purpose 
of identification, the Kalman filter representation is 
desirable 


x 


A A 

Fx + Gu + K(y - Hx) 


( 2 . 12 ) 


y = Hx + v 


[F , G, H , 1 ] have been mapped into [F , G, H , K ] . It is usually 
easier to identify K rather than 1 , and in general 
requires fewer parameters. 


In the state variable form, the natural frequencies and 
damping ratios of all modes are estimated directly. The 
measurements are, however, nonlinear functions of parameters 
leading to a more difficult convergence problem. 


o 



The propagation of the linear equation, as well as the 
computation of the gradients of the innovations with respect 
to parameters are much more easily accomplished with a dis- 
crete space model of the form 


x 


k+1 


= Fx k + Gu k + fw k 


= Hx. + Du, + v, 
k k k 


where w k and v k are independent, identically distributed 
sequences with zero mean. If the time variation of the 
estimator during the transient phase is important or there 
is a priori structural knowledge about the noise covariance matri- 
ces that considerably reduce the number of parameters below pn , 
the dimension the Kalman gain matrix K (where p is the 
number of outputs and n the number of states), then the 
general equations above should be considered. In real time 
flutter identification, the r matrix generally will not be 
known, nor the structure of the noise covariances, making 
the innovation form a much more desirable set of eqautions. 

They are given by 


x 


k+1 


Fx k + Gu k + Kv k 


y k = Hx k + Du k + v k 


where (v k ) is a sequence of independent variables with 
zero mean. Moreover, the model generating y from y , 
i.e., y = II(zI-F+K H) y , is stable (eigenvalues of 
F-KH are in the unit circle). In the next subsection, we 
will see that ARMA type-models have direct realizations as 
a state-space innovations form. 

The discrete controller and real Jordan block form, analogous 

to the continuous block forms described previously, both keep 
a simple parameterization of an oscillatory mode; however, 
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the real Jordan block form is a preferable canonical form 
for identification of physical systems, generally giving a 
better conditioned optimization problem as described below. 

Consider a single mode (a general system is a linear com- 
bination of such modes) of a continuous system 


— # — 



— 


— — 



X 


0 

l 


X 


0 

• 

V 


2 

-0) 

L n 

-2£co 

nj 


V 


1 



y = 

[h i V 




v 


The input is an acceleration forcing the velocity 
equation, with position the integral of velocity. The 
discrete form of this equation (for a pure oscillator using 
the simplest case) is 


V 


1 

*~1 

sin(to n A) 


x k 


r 2 , - 

“n - 1 

J, 'k+1 

= 

X 

w 

n 



4 

% 

v k+l 


sin(« n A) 

1 


v k 


sin( 03 n A ) 


G) 

n 

2 

03 

n 


K 


0) 

n 

. _ 






x. 


y k = [h l V 


y 


k 


If time is normalized such that for the primary mode of 
interest, co = 1 , then the above discrete equations are in 
real Jordan form. Consider, in addition, the identification 
of the mode in both the block real Jordan form and the block 
controller forms of Equation (2.10). The two equations are 
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X 


= Jx, + G u, , x. ... = F 


k+1 ' u m u k ’ ~k+l x c x k ' ~c~k 


+ G u. 


y k H m x k 


y k = H x k 


where in both cases the H matrices have been changed to fix 
T 

G at [01] , and conserve the residues. For our one mode 

example the parameter vectors are 



h 

m l 



’ h l~ 




h m„ 



h 2 


0 1 

9 

z 

and 

0 = 


for 

-b -a 


a 



-b 




Cl) 



-a 



The gradients of the estimate with respect to these two 
different parameterizations expressed in the real Jordan form 
state variables are 


Te = T m !x l’ x 2' x l ( - ) ' x 2 ( - ))T ' 


and 


-^r = T[x 1 ,x 2 ,x 1 (-) ,x 2 (-) ] T , 

8 0 



0 

0 ) 


The Hessian, or information matrix is approximated by the 
outer product of these two gradients. If 
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M 


then 


x 1 (k) 

x 2 (k) 

x 1 (k-l) 

L x 2 (k_1) 


N 

X/ 4* $ 

k=l 


M 


m 


T MT ^ , and M 

mm’ 


= T M T‘ 


In other words the transformations, T and T are 

directly effecting the condition of the estimation problem. 

Since we are referencing both to the real Jordan form, 
f 

K 2 (T m ) = 1 , or it is perfectly conditioned with respect to 
its own coordinate system whereas k 2 (T) can be quite large. 
It was pointed out above that continuous modes driven 
directly by accelerating forces, give discrete models that 
are closer to being in real Jordan form, hence transforma- 
tion to the block controller form can only make the numeri- 
cal conditioning of the information matrix worse. Of course 
in high order physical systems the phase relationships 
between modes can be such that this preference for the real 
Jordan form over the block controller form would no longer 
hold. However, the primary modes in flutter testing are 
expected to sufficiently follow the phase behavior of a 
simple oscillatory mode, that the real Jordan form will be 
used . 

Having described different model forms and reasons for 
preferred ones, the next two subsections describe the 
algorithms used to identify the parameters of these model 
sets . 


t 


The condition of a matrix A with respect to any norm x is 


K x< A > 


Mix IM -1 II x 
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2 . 2 Recursive Identification Algorithms 


Recursive identification refers to identification 
algorithms that update parameters at every sampling instant. 
Such algorithms are characterized by finite non-increasing 
storage requirements. Typically they are well-suited for 
on-line, real-time identification on computers of modest size. 

For single-input single-output systems in linear finite 
dimensional form, the AutoRegressive Moving Average with 
exogenous inputs (ARMAX) models are popular. They are 
represented with the difference equations: 

(l+a lZ - 1 + a2 z- 2 ...+a NA z- NA )y k+1 = (VV” 1 ' • • +b (NB-l) z ~ <N 
+ (l+c 1 z- 1 +c 2 z- 2 ...+o Nc z- NC )e k+1 

or 

A(z _1 )y k+ i = B(z _1 )u k + C(q -1 )e k+1 
where 

y k’ y k+l’ - • * are outputs 


u k ,u k _^,... are the exogenous inputs 


e, , e. are the innovations or white noise, 

k’ k-1’ 

z ^ represents the unit delay operator, i.e., z ^y k = y k 

Inputs and outputs are observed but the innovations 
sequence is not observed. In the transfer function language, 


, -l s -1 B(z 1 ) , -1 N 

y (z ) = z " _y- u(z ) 

A(z X ) 


C(z 1 ) 
ACz- 1 ) 


v(z 1 ) 


C(z ■*") always represents a stable polynomial, i.e., all its 
roots are inside the unit circle. In general B(z and 
A(z are not stable and not comprime. A general black-box 



model 


F(z 1 ) Y( z X ) = 


B(z 1 ) . -1. 0(2 ) 

—~ZT u ( z ) + ZT~ 

A(z ) D(z ) 


can also be used but they have not been very popular, partly 
because the ARMAX formulation can subsume the general form 
without substantial penalty. The ARMAX form can be seen 
as an implementation of the so-called observer cannonical 
state-space form, 



-a 1 10. . . 0 
-a 2 01.... 

~ a 3° \ 

1 

-a 0 0 

n 

H = [1,0. ..0] 

The vector ARMAX models consist of the vector difference 
equat ions 

A (z _1 ) y k+ 1 = B(z _i )u k + C(z -1 )e k+1 

where 


y k’ y k-l 

are 

vector 

outputs , 

u k’ u k-l 

are 

vector 

inputs 

e k’ e k-l • • ' 

are 

vector 

innovations 
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Some special model forms also have names, 


A(z 1 )y 


k+1 


B( Z _1 )u k + e k+1 


is known as equation error form . It is obtained as a special 
case of the ARMAX form by setting C(. •) = 1 . 


Setting P( • ) - 0 and C( • ) - 1 gives the autoregressive 
(AR) form: 

A(z " 1 )y k+i * e k+i • 

Setting A( z b = 1 and E(z b = 0 gives the moving 
average (MA) form: 


k+1 


C(z 1 ) v 


k+1 


The ARMAX forms most useful in identifying flutter modes 
are the recursive least squares (RLS) model (F(z -1 ) = C(z -1 ) 
= D(z -1 ) = 1) , useful with persistent known inputs of 
reasonable magnitude with a very large signal-to-noise ratio, 
and the recursive maximum likelihood (RML) method [2] 

(F(z "S = C(z '*') = 1) , useful for modeling with unknown 
input (the B polynomial can be zero when there is no known 
input) or with considerable noise disturbance. Identifying 
the C(z 1 ) polynomial is equivalent to identifying the 
Kalman gain directly in a state-space innovations form model 


(described previously). 


2.2.1 Recursive Least Squares (RLS ) 

RLS is a very popular recursive identification 
algorithm. It is very simple to implement and provides 
reasonable estimates with excellent robustness. 

RLS gives exactly the same estimates as the batch least- 
squares on the cumulative data. In addition to giving the 
least squares parameter estimate at each time instant, RLS 
also permits forgetting of past data to accommodate quasi- 
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stationary models. It also permits incorporation of prior 
uncertainty about the starting values of the parameter 
estimates. The model underlying RLS is the same as the 
one for batch least-squares, 


y ( k) = (J> ( k ) 0 ( k-1 ) + e(k) 

In the context of linear dynamical systems the auto- 
regressive form and the equation error forms can be trans- 
formed to RLS form with 


(1 + a^z z a n A z = 



1 

1 

'C 

1 

H 1 

. 1 


a i 

II 

X 

-e- 

y (k-2 ) 

e = 

a 2 


~ y (k-nA) 


a NA 


ft . * * 


— - 


( 1+a^z 


. + a 


NA 


-NA s 
z )y 


k+1 


“ < V b l Z_1 ' • • +b NB-l Z '‘ < ™" 1> )u k +e k+l 



Y (k-1) 


a i 


i 

V! 

/ — V 

9? 

1 

CO 


a 2 

il 

X 

-e- 

~ y (k-n) 

CD 

II 

a nA 


U (k-1) 


b 

o 


U ( k-2 ) 


b ! 


U (k-N B ) 


t) 

NB-1 
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2.2.2 Recursive Maximum Likelihood (RML2 [2]) 


Estimation of the A, B, and C polynomials of the 
ARMAX model : 


A(z 1 )y k = B(z X )u k + C(z 1 )e k , 


is called recursive maximum maximum likehood. For SISO 
systems, since A and B are both monic polynomials, this 
equation can be written as 




1 - 


a(z 1 ) 
c(z -1 ) 


■ b(z ) 

y + ' u + q 

, -1. k k 
c(z ) 


The one step predictor of y k based on a certain set of 
values of a^, b^ and c^ , denoted by 9 , is as follows 


.a 


y k|9 


a(z 1 ) \ b(z - 1 ) 

, -lv I y k t -1. 

c(z ) / c(z ) 


U, 


The one-step-ahead prediction error is given by 

e k ( 0 ) - y k - y k l 0 k _! 

The prediction error methods are based on the minimization 
of a quadratic function of the form 

J - §X»> 

k=l 

The above minimization problem is nonlinear in 0 . Hence 
an explicit solution is not available. Therefore, a numeri- 
cal search procedure is used to find 0 , 

Differentiating J with respect to 0 gives 
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N 

J = T] e k (6)Ve k (0) 
k=l 

VJ and Ve k (0) are the gradients of J and e k (0) • The 
second gradient is given as follows 

N 

V 2 J = £ tVe k (0)V e R (0) + V 2 e k (0)e k (0)] 

k=l 


From the prediction error equation 


Ve k (0) = 





2 

The Hessian, V J , will now be approximated at the true 

2 

minimum of J ; the second term of V J can be shown to be 

zero at the true parameter value. Since the true Hession is 

more important close to the true minimum than elsewhere, we 

2 

approximate the Hessian by the first term of V J which, in 
view of Ve (0) , becomes 

rC 

N 

v E Vy k|e v y k | 0 

k=0 


Further approximations are needed in order to obtain a recursive 
Gauss-Newton algorithm from the above equations and compu- 
tation of e j^(0) requires all the data up to t . This 
computation is approximated by using the latest values of 
the data and the parameter estimates. Different approaches 
to affect this approximation lead to several algorithms. 

The RML2 algorithm uses the following approximation 
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e(k) = y(k) - 4>(k) 6 (k-1) , 

where 


<fr(k) 


and e = - cj>^_ 0 ^ are the residuals or the a posteriori 

prediction errors. 

2 . 3 Block Recursive Identification Algorithms 

In block recursive or semi-batch algorithms a block of 
data is stored in active memory and the identification 
routine may make several passes over the data. When the new 
data block is brought in, the old one is discarded. Maximum 
likelihood techniques are well suited to this situation, 
particularly in the known input case where an output error 
minimization can be used. Sensitivities of the outputs with 
respect to parameters are propagated as well as the estimate. 
The parameters of the Kalman gain in the innovations form 
can also be propagated for the unknown input case. 

Two algorithms were developed in a high-level inter- 
pretive language [1] to study their performances. The 
first propagates the analytic sensitivities of a SISO system, 
propagating the minimum number of equations possible. The 

second algorithm is structured for a general MIMO system 
with potential nonlinearities. The sensitivities are found 
by perturbing the nominal model. The following sets of 


y(k-l) 
y(k-n ) 

u(k-l) 

u(k-n b ) 

i(k-l) 

e(k-n ) 
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discrete equations are propagated 


X 


k+l 


f (x k ,u k> 0) 


X 

o 

f(x k ,u k ,e+A8 l ) 

, X 

0 

• 

’ o 


f<x k’ u k, 6 +Ae m ) - 


0 


(m+l)nxl 


where 8=0 


mxl 


The outputs and associated innovations at each time point k 
are : 


z o 


h( x 0 > u , 0) 


"v 1 


z -z 
o 


V 

Z 1 

= 

hCx^u ,0+A0 1 ) 

=> 

V 1 

= 

(z 1 -z)-(z Q -z) 

and 

• 

z m 


h(x ,u, 0+A0 
m m 


_ v m__ 
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v m 


Rather than analytic sensitivity time histories as the 
regressors, the perturbed sensitivities are regressed on the 
nominal innovations to identify a scaled version of the 
parameter step, 

v o = || se = nse/Ae) 

The parameter update equation is 

0£ +1 = - (60/A0)^ • diag(AG^), where l is the 

interation count. The parameters associated with measurement 

equation above are updated with a U-D update performed once 
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for the entire data block. The measurement covariance is 
estimated by 

R £ + 1 = diag{[v o -nse/AO)] [v o -T( 60 /Ae) T ]} . 

Results using this algorithm on a two mode example with very 
lightly damped modes are shown in Section 4.3. 
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SECTION 3 


ALGORITHM RELIABILITY, ACCURACY AND ROBUSTNESS 


The previous section described both a general class of 
recursive estimators called recursive predictive error 
methods and a semi-batch maximum likelihood algorithm for on- 
line flutter identification. Application of such algorithms 
to flutter identification can not be made in a useful way 
without a qualitative and quantitative assessment of the 
reliability, accuracy and robustness of the algorithms. By 
reliability we mean, can the algorithm fail to converge, 
converge to a wrong value or give grossly erroneous accuracy 
estimates? An accurate, but large parameter error estimate 
is reliable information, quantifying the appropriate confi- 
dence in poorly identifiable parameters. The accuracy of 
parameter estimates can be assessed in three ways; simula- 
tion, the Cramer-Rao bound or estimates of this bound from 
either the propagated covariance or limit evaluation. All 
three accuracy measures and their use with flight data are 
discussed in this section. Robustness of flutter identifica- 
tion algorithms addresses the ability of the algorithm to 
maintain reasonable accuracy under violations of assumptions 
on the model and experiment. Thus, all three of these 
algorithm attributes are inter-related as will be elaborated 
in subsequent discussion. 

3 . 1 Algorithm Reliability and Convergence 

The primary reliability issue for flutter identification 
algorithms is convergence. 

• Do the estimates of the modal parameters converge? 
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• If the estimates converge, what is their limit? 

Is it unique? 

• What is the convergence rate and the accuracy of 
the estimates? 

The semi-batch maximum likelihood algorithm requires 
good a priori estimates for convergence, but the convergence 
is quadratic in the neighborhood of the true values. An 
example is shown in Section 4.3 where RPEM applied to only 
100 points can give excellent ML start-up values from an initial 
RPEM estimate of zero for all parameters. Quite a bit can 
be said about the convergence of recursive predictive error 
methods. The conclusions are affected by the user choice 
of model set, choice of input signal, choice of criterion 
function, choice of gain weighting sequence, choice of 
numerical serach direction, and choice of initial conditions, 
as well as the selection of the "approximate" gradient of the 
predictor (see Ljung and Sodestrom [3] for a thorough dis- 
cussion of all these issues). Here we summarize general con- 
clusions on RPEM applicable to the chosen model sets (see 
Section 2.2) and experimental conditions experienced in 
flutter mode monitoring. 

The recursive predictive error method will converge 
with probability one to a local minimum of the expected 
cost as N, the number of observed data, tends to infinity . 
Whether local minima exist apart from the global minimum, 
depends on the model strucure chosen. For the auto regres- 
sive form 

A(z)y k - e k , 

and ARMA form 

A(z)y k = C(z _1 )e k , 

no false local minima exist [4]. For the equation error 
model set , 

A(z _1 )y k = B(z _1 )u k + e k , 
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no false local minima exist since the cost is quandratic in 
the parameters. 

For the ARMAX (Auto Regressive Moving Average with 
exogenous inputs) model set, 

A(z _1 ) y k+1 = B(z -1 ) u k + C(z -1 ) e k+1 

The negative log likelihood function of the ARMAX form has a 
unique global minimum if the signal-to-noise ratio is sufficiently 
large. Local minima do exist if the signal-to-noise ratio 
is very small [5]. This will be the principal model used, 
with the Recursive Prediction Error algorithm (RPEM [3], 
also known as RML2 [2] for the ARMAX form), which implements 
the Gauss-Newton update for parameter estimates. 

The RML1 algorithm [2] and the AML (approximate Maxi- 
mum Likelihood) algorithm of Solo [7] use a further approxi- 
mation in the Gauss-Newton update which in turn poses the 
so-called Strict Positive Real condition for convergence [8]. 
However, the projection of the estimates into a region of 
stability is not needed with AML. Even the problem of the 
estimates getting trapped into local minima disappears. 

If the ident if iability conditions are not met, the parameter 
estimates do not converge, but the prediction-error sequence 
does converge to its minimum value. 

Starting RPEM or RML2 as AML and gradually shifting to 
fully RPEM or RML2 is possible by using the so-called con- 
traction factor on the estimated C(z ) polynomial [3]. 

Doing so eliminates the problem of getting trapped in the 
local minima if they exist. The transient convergence rate 
is also markedly improved. 

Another significant assertion about the above conver- 
gence results, is that convergence holds whether or not the 
actual system has the same model structure as the chosen 
one. The identification procedure will pick the "best" 
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approximation of the system, or in other words the approxi- 
mation from the selected model structure which best mini- 
mizes the prediction errors. However the identified model 
may depend on the input, since the best approximation for a 
sinusoidal input, for example, may be different than for a 
white noise input. If the selected model structure is 
exactly the same as the true system, the "best approximation" 
is equal to the true system, independent of the input. 

In summary, these black-box models are highly desirable for 
identification of flutter modes where a detailed model of the 
aircraft may not be available and could change considerably during 
a flutter test. 

Thus far we have described whether the proposed algorithms 
converge and if so, whether they converge to a unique global mini- 
mum. In the next section a description of asymptotic distribution 
of the estimate 0 ( t ) will be given, which provides both the rate 
of convergence and the accuracy of the estimates asymptotically. 

It should be emphasized that there is no similar analysis for the 
transient convergence behavior of recursive algorithms, and that 
this behavior is very strongly a function of the "user choices" 
mentioned at the outset of this section. Transient convergence 
behavior can best be investigated by simluation. The simulation 
performance, of the model analyzed theoretically in Subsection 3.22, 
for short, intermediate, and long data records is given in Section 4. 

3 . 2 Estimation Accuracy 

The semi-batch maximum likelihood method and the recur- 
sive prediction error technique propagate the covariance 
matrix (inverse of the Cramer-Rao information matrix). This 
matrix is propagated in the U-D factored form. The Cramer- 
Rao bound is based on an assumed noise distribution, typically 
Gaussian noise sources, which may not hold in flight tests. 


26 



Also the bound holds for the asymptotic nature of the 
parameter estimates. The physical environment and system 
characteristics can change before a long enough data record 
can be processed, to achieve the asymptotic behavior. 

The estimation errors predicted by the Cramer-Rao bound 
are usually too small. However, it provides an excellent 
measure of relative parameter estimation errors. Experience 
[9] has shown that in stability and control derivative esti- 
mation, the Cramer-Rao bound should be multiplied by a 
factor of three to five to obtain actual estimation errors. 

In model parameter estimation with high band width sensors, 
the factor is closer to two based on limited experience. 

The simulation cases of Section 4.4 show the actual errors 
and the theoretical bound (computed from the true model). For 
the long data simulation sequences, the bound gives a good 
accuracy estimate for the input noise sequences, which are 
Gaussian in the simulations. 


For both the semi-batch maximum likelihood approach and 
the recursive predictive error method, the information matrix 
inverse is approximated as 


P 



3V £-1 9V T 
3 0 K 3 0 


L -i 

In other words the second gradient of the criterion function 
(V) is approximated as a weighted outer product of the 
first gradient of V with respect to the parameters, where 
R is the estimated innovations covariance. 

3.2.1 Maximum Likelihood Estimation Accuracy 


As explained in Section 2.3 the parameter update in 
the semi-batch ML algorithm which was the inverse of the 
Hession, is accomplished with the Bierman's U-D update 
routines [ 10 ] . 
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The variances of identified parameters can be obtained easily 
from the U-D factored form of the covariance. These parameters 
are, however, parameters of a discrete model used for efficient 
propagation, whereas we are interested in the continuous fre- 
quency and damping estimates and their standard deviations. We 
can use the Jacobian of the transformation from the discrete 
roots (the ML parameterization) to the continuous frequency and 
damping parameters. This is a local linearization of this non- 
linear mapping; however, for lightly damped modes, this trans- 
formation is well-conditioned indicating that it is not sensitive 
to small variations in discrete parameter values. 

Since the parameter covariance is 

P„ = E[A0 A0 T ] , 
u 

the covariance of different parameters which are functions 
of the original ones 

Q ± = f ( 9 ) , 

% " Jp e jT 

where 

J - If * “i - Jie ■ 

It is easier to write the inverse function in this case 
0 = f“ 1 (0 1 ) , 

where 0 are the discrete parameters and 9^ the contin- 
uous frequencies and damping ratios. For each mode, the 
discrete parameters are given by 
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°D = 6 
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-£a)^x 


cos[w n (l - % ) 2 t ] 


sin [ ai ( 1 - ? 2 ) 2 t] 


D - n 

and their variations are given by 
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where x is the sampling time and -co = co (l-£ ) 2 the 


D 


n 


continuous damped frequency. With the U-D factored co- 

variance it is necessary to "square up", the new covariance. 

However, for accuracy prediction, any loss of precision due 

T T 

to "squaring up", is insignificant m P Q = JUDU J , where 

6 1 

for k modes 


r i 


J = 


NH 


parms 


[J L ] 


[J k ] 


I, 


NG 


parms 


This estimate of the Cramer-Rao bound from the inverse of 
the approximated Hessian is accurate only in the region of 
the convergence point . This is demonstrated quite clearly 
at the convergence point on the simulated example of 
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Section 4.3, where the actual estimation errors show nearly 
quadratic behavior as well. 


3.2.2 Recursive Predictive Error Method (RPEM) Accuracy 
Est imates 

The internal covariance estimate propagated in U-D 
factored form in the RPEM algorithm is not suitable to esti- 
mate the Cramer-Rao bound for short data records because 
it is regularized to prevent singularity and also effected 
by the transient effects in the exponential forgetting fac- 
tor which weights recent data more significantly than past 
data. The Cramer-Rao bound can be computed for RPEM algori- 
thms, however, exactly in the simulation case, where the true 
parameters are known, and estimated during a flight test propaga- 
ted UD covariance after the initial convergence transient. 

Consider the ARMAX model, 

A(z _1 )y k = B(z -1 )u k _ 1 + C(z -1 )e k , 

with the negative log likelihood function 



as described in Section 2.2.2. The individual component 


derivatives of 


3y t 


3 *k 

30 

1 


are 


9a i C(z 1 ) Yk 1 


8y 


k = 


u 


9b. C(z - 1 ) 

l 


k-i ’ 
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1 


1 


9y k 

3c. 

l 


Thus gradient of 


C(z _1 ) k ' i 

TT " *k = 


- 1 T k 

C(z X ) k 


The regressors are 


"filtered" versions of the output, input, and innovation 
sequences . 

We seek the asymptotic parameter covariance, where 

N 

H _1 = e[^R"V]= ~ £ E[^R“V] . 

k=l 


and N is the number of data points. 


Just as the ARMAX model above can be realized as a 
state-space innovations model, the gradient can also be 

realized as a linear model driven by the sequences u^ and 
e^ . This realization uses an augmental state and input 



1 

1 


Vi 

11 

X 


w k-l 



1 

-G- 

1 • 


e k-2 ^ 


in the equation, 


F ? k -1 + <3 * k _! 


For the scaler output case, 

P = - , with t = HZH ' , 

NR 

where Z is the solution of the discrete Lyapunov equation 

T T 

Z = FZF + GQG . 
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Q is the covariance of , and H selects the upper cor- 

ner of Z corresponding to the row dimension of , The 
realization of the equation above uses the relations 

(1 + cCz" 1 ))^ = cj) k _ 1 . 

and 

y k = ^k-i^k + e k 

giving the following augmented equation. 



Given quite general conditions on the experiment (see 
Ljung [3, Chapter 4], and conditioned on the event 0^ 0 

with probability 1 (w.p.l), then the constant scaled Lyapunov 
solution 
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NP 


RT 1 (0*) 


is the asymptotic variance of the statistic 



(0-0 )e N ( 0 , NP ) . 


Noreover, for any 6 > 0 

N^ _<S |0(N) - 0*| ->0 w.p.l as N -> oo 


(see Ljung [11] for the general, formal statement and proof 
of this theorem) . 

* 

In simulation cases where 0 is known the Cramer-Rao 
bound above can be computed and used to estimate the theo- 
retically attainable accuracy of the parameter estimates as 
a function of time, whereas with experimental data only the 
identified model can be used to estimate this bound. The 
simulated results of Section 4.4 compare the accuracy of one 
single simulation and the theoretical Cramer-Rao bound for 
different data lengths. 

1'n contrast to the maximum likelihood formulation, 
when the model is already in discrete modal form, an addi- 
tional transformation must be applied, to convert the discrete 
parameter covariance matrix to the continuous frequency and 
damping parameter covariance as follows, 


T T 

P 0 = J X P 0 X J , 


with J defined as in the previous subsection, and 
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M is the modal matrix which transforms the polynomial 
parameters 



The above transformation was applied to the Cramer-Rao bound 
based on the true simulation parameters to produce the results of 
Tables 4-3 through 4-5. A high eigenvalue condition number, i.e., 
s^(A_^) = ^ X i y i® (norm of A^'s spectral projector; x^, are 

the left and right eigenvectors), is a warning that the respective 
polynomial roots are very sensitive to errors in the polynomial 
coefficients. For the simulation results described in the tables 
mentioned above, the eigensystem condition for both the true 
parameters and the estimated ones was a reasonable level, less 
than approximately 300 in the worst case. 

In the known input case, which is the recursive least 
squares, also called the equation error form, 

A ( z_1 ) y k = B(z_1)u k-1 + e k * 

the gradient is 



80 

where 

^k = [y k-l y k-na ’ u k-l •** U k-nb ] 

For arbitrary inputs u^ , the asymptotic covariance is 
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-1 T 

V ] 


1 


This parameter covariance can then be found in terms 
of autocorrelations and cross-correlations of the input and 
output sequences. 


For the special case where u and e. are random 

k k 

sequences we can compute the parameter covariance directly 
by solving a discrete Lyapunov equation as in the ARMAX 
case, which employs the covariance of the two random 
sequences. The state space realization uses the state and 
input , 



in the equation 


K 


k 


F? 


k-1 


+ Gw 


k-1 


with F and G given by 



Computation of the Cramer-Rao bound as described in 
this subsection allows the investigation of the effect of differ- 
ent gust and sensor noise levels without performing multiple 
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simulations. This is discussed in the next section under the 
general topic of robustness. 

3 . 3 Robust Flutter Mode Identification 

Violations in the assumptions of the statistical nature 
of the disturbances require estimation techniques robust 
with respect to distribution. These techniques are summa- 
rized and demonstrated by simulation. Violations in the 
assumptions of the number of active modes, the gust levels, 
or characteristics of the sensor noise effect the accuracy 
of estimation. The results of the previous section can be 
used to efficiently quantify these effects. 

3.3.1 Robustness with Respect to Distribution 


Robust estimation procedures are very important 
because a few bad data points can significantly degrade 
parameter estimation accuracy. In some cases convergence 
characteristics would also be affected. Ljung [3] summa- 
rizes robust estimation problem as follows: 

1. When the measured data set contains some values 
that are abnormal, e.g., quite large due to 
sensor failures, straightforward use of a quad- 
ratic criterion function will give substantial 
jumps of the parameter estimates. Moreover, 
the estimates need a long period before con- 
verging back to the previous levels. 

2. A way to cope with this problem, i.e., to 
robustify the algorithm, is to use a criterion 
function that grows more slowly with e than 
the quadratic one. Then large prediction errors 
will get less influence on the parameter esti- 
mates and the algorithm becomes more robust. 

3. Another approach is to test recursively if the 
data contains outliers. This can be done by 
comparing the prediction errors with a speci- 
fied limit. Large prediction errors mean that 
an outlier or a measurement error is probable. 
The measurement can then be substituted with the 
predicted value. This approach is applicable 
when there is only a few outliers in the data. 
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Both identification techniques we have discussed in this 
section are based on minimizing a quadratic penalty func- 
tional, e.g. 


N 



k=l 


The robust procedures modify this cost functional to include 
a weighting function g which is monotonically decreasing 


V 


sL 


k=l 


e 2 

k 2 

B k 
k 


where B, is the covariance of e, 
k k 

The use of such a robustness procedure in the real-time 
flutter identification procedure will now be discussed. 


Recursive identification techniques are inherently 
well-suited for implementing robust estimation techniques 
because the reasonableness of each new data point can be 
evaluated before it is used to update the parameter esti- 
mates. An illustrative discussion of robustness measures 
and robust identification is given by Ljung [3]. A sum- 
mary of the steps necessary to make the parameter estimates 
robust against data dropouts, wild points and outliers is as 
follows : 

Step 1: After a data point is processed, the 

estimated parameter values are used to 
predict the next data and the standard 
deviation of the predicted value. Let 
the predicted value be y with stand- 
ard deviation a . ^ 

Step 2 : The next data point is compared to its 

predicted value. The quantity 
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£ = 


y - y 


p 


is a measure of the probability that 
the data point is wild. Data points 
with £ > 3 should be suspect. Thus, 
data corresponding to large values of 
£ should be given less weighting. 
Depending on the expected distribution 
of data droupouts, any of the weighting 
functions of Figure 3-1 can be selected 
(see Huber [12]). if f(£) = y^ , the 
estimates have the robustness of a median. 


Step 3: A data point completely outside the predicted 

distribution is not used 

A restart is necessary if several contiguous points are re- 
j ected . 



£ 


Figure 3-1 Weighting Functions for Enhanced Robustness 

A weighting function of the second type was demonstrated, 
where the cost function is 


e <a 
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The measurements with predicted errors out of bounds are replaced with 
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Data dropouts have been generated with distribution shown 

below . 


Amplitude 

Normally distributed 
variance as the data 
at 5o from the output 

with the same 
but centered 
est imate 

Frequency 
and sign 

Semirandom telegraph 
X = 1% 

signal with 


Results with and without these dropouts are compared in Table 3-1 
for the known input case. (This table with a 6-model identified 
mode can also be compated to the same simulation in Tables 4-3 
and 4-4, which have no data dropouts and use an 8-mode model.) 


o> n V ~ "n U%) ^ (%) 


125 points 
with robust 

1 

.9998 

4.70 

8.50 

est imat ion 

2 

1.9592 

4.80 

6.36 

125 point 
w/o robust 

1 

.9998 

4 ) 70 

8.48 

estimation 

2 

1.9597 

4.80 

6.39 

500 points 
with robust 

1 

.9997 

4.70 

4.78 

est imat ion 

2 

2.0009 

4.80 

4.71 

500 points 
w/o robust 

1 

.9980 

4.70 

4.87 

estimation 

2 

2.0037 

4.80 

4.51 


Table 3-1. RPEM Frequency /Damping Estimates ( With and Without 

Robustness Techniques - Known Input with Standard 

Gust Cases o /a gust = 25) 
u 

The robust estimation procedure can improve the accuracy of 
the flutter parameters (see Table 3-1); however, during the 
rapid portion of the algorithm convergence, the wild points 
actually increase the convergence rate. The extent to which the 
model is overparameterized also has a significant effect in 
making RPEM robust. 8-mode models identified from this data 
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showed little degradation until asympotic accuracies were 
achieved, and then the degradation was not severe. 

3.3.2 Effects of Noise and Modeling Errors 

Very high frequency effects can be modeled as white 
noise, unmodeled system dynamics with colored noise, and 
low frequency effects as a random walk. Errors in estimation 
are caused by an incorrect order linear approximation model, 
an incorrect intensity of Gaussian noise, or finally assump- 
tion that the noise distribution is, in fact, Gaussian. 

The effects of these different noise sources on the 
Cramer-Rao bound for the modal parameters can be computed 
very efficiently for a large number of cases by merely solv- 
ing matrix Ricatti and Lyapunov equations for each desired 


situation, rather 

than performing 

a simulation and doing 

the identification for a long data record. 

Consider simu 

lations generated 

for the 

2 -mode 

model of 

Section 4.4. 

according to 

the 

following 

model , 



X k + 1 

= 

Fx k + Gw k 

+ Gu k 



y 

= 

Hx k + V k 

9 



with 






W k 

= 

N( 0 , q ) , 

q * 

2 

°gust 


U k 

= 

N(0,ajj) 




V k 

= 

N( 0 , r ) , 

r = 

2 

a 

sensor 

f 


but the sequence u^ is known. u k can be considered 
wideband excitation which could quite possibly be a pseudo 
random binary sequence; we have made it Gaussian here for 
computational simplicity. 
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The ''nominal” model, about which trends for how these 
different disturbances or inputs effect the achievable 
parameter estimation accuracy will now be described. The 
two modes have frequencies of 1 and 2 per rev (borrowing 
a term from rotor craft analysis for normalized frequencies), 
where 1 per rev is 11.8 hertz-emulating the open loop modes 
in the ARW-1 vehicle. The sampling rate is approximately 300 
hertz. The gust noise covariance (a - " . = (.15) deg/rev^) 

gU S L 

has been chosen to correspond to a continuous model with vertical 

4 * 

gust root-mean-square values of approximately 1-2/3 feet/sec (at 
the ARW-1 condition of Mach of .7, altitude of 15,000 feet, and a 
correlation distance of 1750 feet). The damping parameters were 
chosen to be similar to the ARW-1 flight condition analyzed in 
Section 5.1 (an 11.8 Hz mode with 4.70% of damping and a 23.6 Hz 
mode with 4.80% damping). The sensor noise observed on the 
data (see Figure 5-11) was taken to be approximately R q = .0050. 
It should be pointed out that the significant difference between 
this simulation and the flight maneuver is the use of broad 
band rather than narrow band excitation. 


The theoretically achievable parameter accuracy as a 
function of input and noise power is shown in Figures 3-2 and 
3-3 for damping ratio and natural frequency estimation errors. 

The flutter parameter variances are a function of the gust, con- 
trol input, and sensor noise covariances, variations of which 
are shown with families of curves in the upper two plots of the 
figures. The three-dimensional plots show the flutter parameter 
standard deviations as a function of power ratios (with a log 
scale on all axes). The parameter covariances have been computed 
via a discrete Lyapunov equation as described in subsection 3.2.2 
they require the control and innovation covariances as input 
parameters. The innovation covariances were computed directly 
from the gust and noise covariances via the discrete algebraic 
Riccati equation [13]. 
t 

'The gust level was very low during the ARW-1, Flight 3. 
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Figure 3 2. Effects of Gust Noise, Sensor 
Input Power on Damping Estimation Errors 
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Figure 3-3. Effects of Gust Noise, Sensor. Noise, and Broadband 
Input Power on Natural Frequency Estimation Errors - Two Mode 

Model 


Note: 3-dimensional axes are Log^Q of the data. 

x-axis is ° u 2 /° w 2 > y-axis is Q w /R ; z-axis is a damping ra tio. 
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The variations about nominal noise covariances were 


' 9 

1 

, and 

R 

4 

1 

1/9 


o 

1/400 


The gust variations. correspond to a range of gust velocity of 
h 2 

Tjft/sec - 5ft/sec , the nominal at Wft/sec for the low gust 

y o 

values seen in the flight data, the low corresponding to negli- 
gible gust, and the 5ft/sec the Dryden gust model. The sensor 
nominal was varied to twice the noise standard deviation and twenty 
time less (corresponding to inertial grade instruments). The in- 
put power variations were 


u 


0 

1 

156.25 

625 


q 


w 


o 


The highest level gave an amplitude response similar to the high-level 
level flight data inputs, and the next lower level is one-half 
the amplitude corresponding to the low-level flight data inputs. 

These theoretical predictions are compared to simulation 
results in Section 4.4. 
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SECTION 4 


SIMULATED DATA RESULTS 

Three groups of flutter examples were simulated to test the 
performance of the RPEM and the maximum likelihood (ML) algo- 
rithms : 

1. A two-mode example with very lightly damped modes at 
a frequency ratio of two. 

2. Different combinations of two modes with light to 
moderate damping ( 1% - 15%). 

3. Two-mode examples with different noise and modal dis- 
turbances discussed in the previous section on robust- 
ness, reliability and accuracy. 

4.1 VERY LIGHTLY DAMPED MODES 

The linear model used to test both the RPEM algorithm as 
the start up estimator, as well as the ML algorithm is 

x = Fx + Gu 

y = Hx 

where 


II 

o 

1 

0 

0 

G = r 0 

-1 

-.0136 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

-4 

-.01 

, L 1 - 


H = [-1.9456 1.2239 2.1867 -.07552] . 


This model was configured to emulate the torsion response (wing 
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bending and torsion modes) of a wing with two modes close to the 
flutter condition with time measured in normalized units (revs). 
The frequency response of this model is shown below in Figure 4-1. 



FREQUENCY 



FREQUENCY 


Figure 4-1. Two Mode Lightly-Damped Model Frequency Response. 
(Frequency in Per Rev) 
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As mentioned in Section 2, it is much more efficient to 
propagate linear systems in discrete form. The discrete form of 

this model (for a sampling time of .23 revolutions) is given by 

/ 


X k + 1 = Vk + Vk 


y 


k 


H D X k 


where 



'. 9723 

.2268 

0 

0 


"o” 


.2268 

.9723 

0 

0 


1 

f d = 

0 

0 

. 8958 

.4420 

’ g d 

0 


0 

0 

-.4420 

.8958 


1 

h d = 

[-.4733 

.2268 

.5027 

-.1105] . 



The 

gaussian 


output of this linear system as forced by a known 
input gust sequence is shown in Figure 4-2. 


47 



st Accel. (Deg/Rev ) Tor. Accel. (Deg/Rev 



Figure 4-2. Lightly Damped Model Output and Input. 
(a gust = * 15 D eg/Rev 2 ) 
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4.2 RPEM AS A STARTUP ESTIMATOR 


One of the desirable characteristics of the RPEM algorithm is 
the ability to start the estimation with no a priori information 
on the parameters. On short to intermediate data records with 
both process and measurement noise, however, overparameterization 
of the model aids in convergence of the flutter parameters con- 
siderably (see Section 4.4, where an 8-mode model was used to get 
excellent flutter identification with two active modes) . In this 
section, where the simulation has no noise, the startup model is 
so good that overparameterization cannot measurably improve it. 

See Table 4-1 and Figure 4-3. Note that the residues of the art- 
ificial mode are negligibly small. 


Table 4-1. RPEM Model Residue Comparison 


True Residues 

4-Mode Model 
Identified Residues 

6-Mode Model 
Identified Residues 

2 . 2683D-01 
-4 . 7332D-01 
-1 . 1049D-01 
5 . 0268D-01 

2 . 2683D-01 
-4 . 7333D-01 
-1 . 1049D-01 
5 . 0270D-01 

-6 . 566QD-09 
-1 . 8139D-09 
2 . 2683D-01 
-4 . 7332D-01 
-1 . 1049D-01 
5 . 0268D-01 
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and y (6-mode model) y and y (4-mode model) 



Figure 4-3. RPEM Performance Without and With 
Overparameterization . 
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These residues of the discrete transfer function correspond 
to the cosine and sine coefficients in a partial fraction expan- 
sion of a continuous transfer function. In other words, for 

_i k 

y(z = H(zl-F) G = E 

i=l 


R c (z+a.) + R g a). 

i l 

(z+a ^) 2 + an 2 


c . and R 
l s . 

l 


are the residues associated with the i-th mode where 


H 


[10 ... 0 ], 


! 1 



. » 
-a , 

G = 

b 

i ; i 



1 0 

y 


. i _ 


- 


and 

H(zl-F) -1 G = H (zI-J) -1 G 
m m 


where H , 
m’ 

J, 

and 

G 

m 

are 

the modal form of the discrete 

equat ions . 

For 

the 

i-th 

mode 

the residues are 


R 

c . = 

l 

h li 

g il 

+ h 2i 

g i2 


R 

s . 

h li 

g i2 

h 2i 

g il 


i 


4.3 MAXIMUM LIKELIHOOD ESTIMATION WITH RPEM STARTUP PARAMETERS 

To demonstrate the use of RPEM as a startup algorithm for 
maximum likelihood the first 100 data points were processed with 
a 3-mode (overparameterized) RPEM model. These parameters and 
the estimated output are equally as good as the outputs shown 
in Figure 4-4. These RPEM parameters are too accurate to reason- 
able test the maximum likelihood algorithm; therefore, an RPEM- 
identified set of parameters using the wrong input (shown in 
Table 4-1) were actually used. 
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As would be expected from the results of the previous sec- 
tion, the fit error of the final model is small, but the parameter 
estimates have not yet converged. 


3 . 


y 6. 

3. 


0 . 


-3. 

- 6 . 

-9. 



0. S0. 40. 60 . 80. 100. 

T I ht! (NORMAL I ZED) 


Figure 4-4. RPEM Start-Up Model. (Estimated from the first 

100 points, but using an incorrect random input.) 


e H 

- 0.4733 

0.2268 

0.5027 

- 0.1105 


- 1.0679 

- 0.7230 

0.3292 

0.6849 


a i 

w i 

°2 

w 2 

0p 

0.9723 

0.2268 

0.8958 

0.4420 


0.9737 

0.2286 

0.8985 

0.4411 

i 


Table 4-2. Maximum Likelihood Start-Up Parameters. 

The discrete 0- parameters have reasonable startup values, 

r 

however the mode shapes are essentially unknown. 
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The; batch or semi-batch maximum likelihood algorithm de- 
scribed in Section 2.3 was applied to this 500 point simulation 
using the start-up values in Table 4-2. Figure 4-4 shows, first, 
the relative error in the natural frequencies on a linear scale 
and then the same error moduli are plotted on a logarithmic scale. 
The estimated variances of these parameters are also plotted with 
the actual errors. The algorithm has converged in four iterations. 
This particular simulation does not have any measurement noise 
added, which leads to the underestimate of variance after the 
modified Newton Raphson optimization has converged. 


53 



Natural Frequency Natural Frequency 

Error Modulus Relative Error 



1 . cl. 3 . 4 . 5 , 6 . 7 . 

Iteration Number * 



Iteration Number * 


Figure 4-4. Maximum Likelihood Natural Frequency 
Errors (and Error Estimates) 


* Note that the last iteration started from iteration five 
with the covariance reinitialized to large uncorrelated values. 
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The damping estimates show similar behavior to the fre- 
quencies, converging in four iterations, with reasonable error 
estimates at the convergence point. Further iterations decrease 
the parameter errors, but without measurement noise on the 
simulation, the error estimates underestimate the parameter error 
on iterations beyond the fourth iteration. 
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Damping Ratio Error Modulus Estimated & True Dumping Ratios 




Iteration Number * 


Figure 4-5. Maximum Likelihood Damping Ratio Errors 
(and Error Estimates) 


* Note that the last iteration started from iteration five 
with the covariance reinitialized to large uncorrelated values. 
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Tor. Accel (Deg/Rev 


Figure 4-6 shows the Torsion rate measurement and estimated 
output with the associated innovations or errors produced by the 
estimated model after six iterations with the covariance correction 
(shown as the seventh iteration on Figures 4-4 and 4-5). 





4.4 RPEM IDENTIFICATION FROM SHORT, INTERMEDIATE AND LONG DATA 
RECORDS 


To demonstrate the convergence, convergence rate and accuracy 
predictions of Section 3.3, simulations with two different models 
were performed. The input was treated as both known and unknown. 
Frequencies of 1 and 2 per rev (normalized frequency) were simu- 
lated in a model of the form 


0 1 



0 

2 or 2 

— oj ^ 2E, ^ 


x + 

1 

0 

1 


0 

2 

W 2 

or 2 

2£ 2 w 2_ 


1 

L J 


where 


y = [0 1 0 -.5] X. The frequency response for the model 
(uj =1, = .0470, oj 2 = 2, E, 2 = .0480) is shown in Figure 4-7. 

This is the same model used to study noise effects on achievable 
accuracy in Section 3.3. The same measurement white noise 
covariance (see Section 3.3) was used (a v ~ .0707), and 
gust covariance (c gust = .05 deg) was used. 

Modal parameters were identified with RPEM as a stand alone 
algorithm for short, intermediate, and long data records, rather 
than merely as a startup algorithm for semi-batch maximum likeli- 
hood. The overparameterized model had eight modes; the lowest 
frequencies are the simulated modes. Frequency and damping 
estimates and the Cramer-Rao bounds (computed from the true 
model) are compared to actual errors in Tables 4-3 through 4-5. 
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Figure 4-7. 


Two Mode Simulation - Damping ratios of 4-7 0 % (for w 
and 4.80% (for w=2). 


1 ) and 






CO 

n 

/\ 

OJ 

n 

°co 

n CR 

?(%) 

t(%) 


Known 

Input 

1 

.962 

7.66lo 

4.70 

10.4 

.0784 

2 

1.896 

3 . 56^0 

4.80 

.94 

.174 

Unknown 

Input 

1 

1.139 

.124 

4.70 

45.2 

12.8 

2 



2.226 

— 

4.80 

93.7 

23.6 


Table 4-3. RPEM Frequency /Damping Estimates and 
Accuracy Estimates 

(2-mode simulation of 125 data points. 8-mode model.) 



0) 

n 

/\ 

0) 

n 

c co 

n CR 

€(%) 

l(%) 

Or % 
^CR 

Known 

Input 

1 

.9986 

3.83lo 

4.70 

4.34 

.0392 

2 

2.0101 

l*78io 

4.80 

5.32 

.0868 

Unknown 

Input 

1 

1.1611 

.0622 

4.70 

31.8 

6.39 


3.525 

.232 

4.80 

11.9 

11.82 


Table 4-4. RPEM Frequency /Damping Estimates and 
Accuracy Estimates 

(2-mode simulation of 500 data points. 8-mode model.) 
































w 

n 

to - 0 ) 

n n 

n CR 

5(%) 

l(%) 

°F °/ 
XR 

Known 

Input 

1 

.99963 

1.7173 

4.70 

4.7083 

.0175 

2 

1.9990 

7 .95 

4.80 

4.6663 

.0388 

Unknown 

Input 

1 

.9798 

.0278 

4.70 

9.35 

2.86 

2 

3.366 

.104 

4.80 

7.90 

5.26 


Table 4-5. RPEM Frequency /Damping Estimates and 
Accuracy Estimates 

(2-mode simulation of 2500 data points. 8-mode model.) 


The conclusions from these simulations both validate the 
performance of the algorithm RPEM for flutter mode identification 
with multiple modes and give insight into the validity of the 
results identified from the flight data in the next section. (The 
flight data also has two active modes with similar characteristics 
and also was overparameterized with an 8-mode model.) 

The gust level is too low to identify the flutter modes with 
turbulence excitation anly. With a known input the algorithm 
finishes the convergence transient in about 750 data points or 
1.5 seconds (see flight data results in Figure 5.15). For long 
data records (>4-5 seconds) RPEM gives extremely accurate esti- 
mates (|£ - £ | < .05%). For this particular simulation the damping 
estimate error (w=l) is less than the Cramer-Rao bound, and for 
w==2, the damping estimate error is less than 2d™ 
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SECTION 5 


FLIGHT DATA RESULTS 


Flutter-test flight data from the third flight of the ARW-1 
vehicle of NASA's drones for aerodynamic and structural testing 
(I)AST) program has been used to demonstrate the real-time flutter 
estimation algorithms. Three different maneuver types were 
processed : 


1. A frequency sweep maneuver with significant excitation 
amplitude , 

2. An intervening record between excitation inputs to 
emulate typical unknown input or turbulence excitation, 
and 

3. The sequence of the last three pulses before the flutter 
incident resulting in failure of the right wing. 


5.1 KNOWN INPUTS 


Symmetric excitation with a logarithmic sweep from 10 -- 40 
hertz is shown in Figure 5-1. 
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The actual aileron positions include some asymmetric input 
due to the Flutter Suppression System (FSS), as can be seen in 
Figure 5-2. 



25588. 25530. 25592. 25534. 25536. 25538. 25600. 

TIMECSEC) 



25588. 25590. 25592. 25594. 25596. 2S598. 25600. 

TIME(SEC) 


Figure 5-2. Symmetric and Asymmetric Inputs 
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The frequency responses of this input and the ALFSO accelerometer 
output are shown in Figure 5-3. The ratio of these two responses 
in Figure 5-3 is a coarse estimate of the transfer function, 
shown in Figure 5-4. 




As can be seen from Figure 5-4, the noise is so great without 
averaging, that there is little recognizable information in any 
frequency range. 

A single-input-single-output identification study was per- 
formed using the symmetric input and the ALFSO accelerometer with 
RPEM. Model orders of three, five and eight modes, with and 
without identifying the C(z _1 ) polynomial (or equivalently the 
Kalman gain) are compared in Tables 5-1 and 5-2. 


Table 5-1. Known Input Modal Estimates 

(Estimation with C(z Polynomial) 


Mode 

3 

Mode Model 

5 Mode 

Model 

8 Mode 

Model 


OJ 

n 

(Hz) 

S(%) 

o n (Hz) 

5(%) 

a) ( Hz ) 
n 

?(%) 

1 

ii 

. 70 

5.77 

11.81 

4.59 

11.82 

4.62 

2 

66 

.72 

22 . 7 

57.19 

11.8 

48.63 

4.00 

3 

228.9 

6.00 

116.8 

17.8 

81.94 

49.0 

4 




173.6 

15.5 

88.53 

11.6 

5 




207.1 

5.42 

143.4 

9.59 

6 






163.8 

5.15 

7 






206.2 

6.28 

8 






238.9 

3.00 



0 . 30 . 60 . 80 . 180 . 150 . 180 . £ 10 . £ 40 . £ 70 . 

HERTZ 


Figure 5-4. Coarse Estimate of A ( j on ) / <5 (joo) 

Z cl 
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Mode 

3 Mode 

Model 

5 Mode Model 

8 Mode Model 


a) n (Hz) 

5(%) 

, 

o n (Hz) 

?(%) 

o n (Hz) 

?(%) 

1 

12.13 

11.1 

11.89 

4.95 

11.77 

4.70 

2 

73.96 

38.2 

52.40 

5.87 

47.59 

5.71 

3 

205.4 

13.1 

113.5 

16.0 

65.44 

18.0 

4 



172.4 

9.24 

103.9 

18.1 

5 



227.2 

8.31 

130.9 

15.1 

6 





16 3.7 

6.94 

7 





198.3 

4.92 

8 





234.3 

3.95 


Table 5-2. Known Input Modal Estimates 

(Estimation without C(z~; Polynomial) 


Only the lowest mode is within the range predicted by 
NASTRAN analysis and ground vibration test. The frequency and 
damping estimate of this mode is comparable and reasonable for 
all six model structures identified except for model with the 
smallest number of parameters, the three mode model with no 
C(z ) polynomial. 

The mode at approximately 11.8 hertz and damping of 4.6% 
is close to the predicted value of the symmetric first body 
bending mode at this flight condition (M = .7, h = 15 k feet). 
Gilyard and Edwards [14, Fig. 3] show predicted values of approx- 
imately a) = 11.8Hz with a damping ratio of 4.7%. 

Note that here we are using the symmetric input 

( <S + )/2 as the single input and hence identifying the open 

L a R 

loop plant even though there is some asymmetric input. The 
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asymmetric input power is approximately 28 db down from the 
symmetric input (see Figure 5-2). Models with the C(z 1 ) poly- 
nomials provide a more desirable model structure, able to model 
this unknown input . 

There are several conclusions on the use of RPEM for real- 
time flutter analysis with known inputs: 

1. Overparameterization is important to accurately 
identify the active flutter mode. 

2. Overparameterization is even more important if the C(z ) 
polynomial is not identified with a known input. 


5.2 UNKNOWN INPUTS 

One desirable characteristic of a real-time flutter identi- 
fication system is to identify the dominant system dynamics 
without a known input for continuous test monitoring. This can 

be accomplished with the RPEM algorithm by identifying the 

-1 -1 -1 
A(z ) and C(z ) polynomials and letting B(z )=0, or equiva- 
lently identifying the system dynamics matrix and the Kalman gain 
directly. Figure 5-6 shows a section of data early in the flight, 
between a sweep and an excitation pulse. 
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The ALFSO sensor output for this interpulse section is 
shown in Figure 5-7. 


855537 . 


85598 . 


85593 . 

TIKE 


85600 . 



Figure 5-7. Turbulence Excited Accelerometer Output 
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Two thousand data points were used to estimate the 48 
parameters and modal estimates of Table 5-3. 


Mode No. 

o) n ( Hertz) 

s(%) 

1 

26.08 

70.1 

2 

50.07 

4.38 

3 

68.57 

20.3 

4 

140.3 

4.31 

5 

171.1 

2.71 

6 

207.8 

8.13 

7 

235.8 

1.31 


* 


Table 5-3. Unknown Input Model Estimates 

(ALFSO Accelerometer Used as Recorded) 


Besides the peaks at 26 and 50 hertz identified above, the ALFSO 
frequency response (see Figure 5-8) shows peaks at 12, 15, 20, 
and 42 hertz as well. 


* Eighth discrete mode not realizable as a real continuous 
system. 
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FREQUENCY (HZ) 

Figure 5-8. ALFSO Accelerometer Spectrum Forced by 
Turbulence 


This estimation was affected by the bias in the acceler- 
ometer, which when taken out gives the modal estimates of Table 
5-4. 


Mode No. 

o n (Hertz) 

?(%) 

1 

14.51 

24.0 

2 

46 . 39 

13.6 

3 

51.43 

7.92 

4 

87.79 

8. 99 

5 

105.5 

4.13 

6 

149.1 

6.94 

7 

191.3 

3.94 

8 

233.4 

.46 


Table 5-4. Unknown Input Model Estimates 

(ALFSO Bias Removed Without Estimation 
of B(z~^) Polynomial) 


The lowest frequency is reasonable, but its damping is too high. 
The estimate errors and error spectrum are given in Figure 5-9. 
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Figure 5-9. Estimate Errors and Error Spectrum 

(A(z -1 ) and C(z Polynomials Estimated) 
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It is evident that the effect of the low frequency modes has not 
been completely removed. The approximate signal-to-noise ratio is 

SNR = (a 2 - a 2 ) / c 2 = 5.77 , 
y e e ’ 

so that the turbulence level is quite low. 

It is interesting to examine the excitation channel more 
closely (see Figure 5-10), which shows what could have been noise 
in digitization, or noise actually transmitted and passed through 
the control system. 
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Figure 5-10. Interpulse Excitation and Excitation Spectrum 

The A, B, and C polynomials were identified using this 
input (a 48 parameter model) and the resulting modal parameters 
are given in Table 5-5. 
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Mode No. 

o (Hertz) 
n 

S(%) 

1 

13.02 

9.39 

2 

48.81 

5.90 

3 

65.33 

16.0 

4 

103.1 

15.1 

5 

137.6 

10.2 

6 

180.0 

7.09 

7 

221.0 

6.81 


* 


Table 5-5. Unknown Input Modal Estimates 

(ALFSO Bias Removed, with Estimation 
of B(z 1 ) Polynomial) 


The lowest natural frequency estimate is still reasonable for 
this flight condition (M = .7 , h = 15 k feet) with the damping 
estimate much more reasonable, for the first bending mode. 

The spectrum of the estimated model errors (Figure 5-11) 
does look like nearly white noise, removing the peak in the 
region of the first symmetric mode. 

The estimated signal-to-noise ratio is slightly larger at 
SNR » 5.79. 


* Eighth mode not realizable as a real system. 
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Figure 5-11. Estimate Errors and Error Spectrum 

(A(z _1 ), B(z _1 ), and Polynomials 

Estimated) 
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The conclusion from this flight data with a very low level 
of turbulence or input excitation is consistent with the con- 
clusions on simulated data in Sections 4.4 and 3.3.2. The ex- 
citation level has a very significant effect on estimation accu- 
racy and, of course, the proportion of the input signals that are 
known significantly affect the parameter accuracy as well. 

5.3 FLUTTER INCIDENT RESULTS 

During the third flight of the ARW-1 program, the Flutter 
Suppression System (FSS), the on-board stabilizing control sys- 
tem, was actually operating at one-half the nominal gain. During 
the acceleration from M = .80 to M = .825 the first wing- 
bending asymmetric mode showed increasingly lighter damping, 
finally resulting in structural failure of the right wing, and 
together with a partial failure of the parachute recovery system, 
a crash of the vehicle [14]. 

The excitation (FSSEXC channel) and ALFSO accelerometer 
response are shown in Figure 5-12. The last three excitation 
pulses are full cycle sine wave doublets (chopped here because 
only every 20-th point has been plotted). Just over twelve 
seconds (6290 data points) were processed as a single record, 
up to the last point before accelerometer saturation. 

Figure 5-13 shows the aileron positions for this final part 
of the flight, while Figure 5-14 shows their average sum and 
difference for the remaining linear portion of the flight. 

It is noticeable that the symmetric and asymmetric inputs 
have comparable input magnitudes due to the flutter suppression 
asymmetric feedback attempting to stabilize the anti-symmetric 
modes . 

At approximately T = 26138.5 seconds, the test pilot was 
instructed to terminate the test and the throttle was reduced 
clue to the light damping observed in the previous time traces 
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Figure 5-12. Flutter Incident Excitation and Response 
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Figure 5-13. Aileron Positions for the End of the 

ARW-1 Flight 3 
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Figure 5-14. Inputs for the Last Linear Portion 
of the DAST Flight 3 
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[14]. The real-time damping estimator indicates a significant 
reduction trend in damping, approximately two seconds before 
this, as seen in Figure 5-15. Such a serious warning, this 
early, is quite difficult if not indiscernable from the time 
trace shown in Figure 5-12.. The last damping estimate (modal 
parameters were output every 100 points or .2 seconds, using 
the last 190 points for the last estimate) is £ = -.02218 for 
the symmetric wing bending mode at 19.16 Herz. Before the drop 
in frequency at the flutter condition (see Figure 5-15), the 
flutter mode frequency was estimated as w n = 125 rad/sec, which 
correlates exactly with post flight processing results by other 
techniques (see Gilyard [14, Figure 12]). 

Figure 5-15 shows the initial convergence period for RPEM 
of 0.5 - 1.0 seconds. It is possible that the damping dip at 
T - 26132 seconds is still due to convergence since some of the 
other discrete modes modeled (an eight mode model) changed at 
this point (after 1500 data points). The break up in natural 
frequency and clear reduction trend in the damping at T = 26136.5 
seconds is definitely due to the changing flutter characteristics. 
All the modes modeled in the overparameterized system are changing 
smoothly at this point. 

It should be pointed out that the asymptotic forgetting 
factor in RPEM can be fixed at less than one to make the algo- 
rithm track models with changing characteristics with higher 
accuracy. Giving the operator the option to select a past 
window size of considered data is a possibility when real-time 
implementation is performed. 
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Figure 5-15. Estimated Frequency and Damping of the Flutter Mode 

(End of the ARW-1 Flight 3) 








SECTION 6 


NEAR REAL-TIME FLUTTER ANALYSIS (NRTFA ) 
PROGRAM DOCUMENTATION 


The delivered ANSI FORTRAN-77, Near Real-Time Flutter Analy- 
sis Program (NRTFA) is so named because, while providing a working 
tool to analyze real time flutter analysis algorithms, it lacks 
the real-time executive, real-time input/output buffered inter- 
faces and real-time graphics drivers. Thus the program reads a 
large data block, then emulates real-time operation by repeatedly 
outputting new frequency and damping estimates to the operator 
console screen. This section 

1) Explains the choice of identification algorithm imple- 
mented in NRTFA,’ 

2) Documents the program structure and subroutines, as 
well as fixed algorithm initialization parameters, 

3) Describes the memory resources and timing comparisons 
for typical flutter models, 

4) Discusses extensions necessary for actual real-time 
implementation, and 

5) Gives installation instructions. 

User information is documented in Appendix A, and Appendix B 
gives the FORTRAN listings. 

6.1 THE NRTFA ALGORITHM 

Studies documented in the previous three sections support 
the conclusion that the Recursive Predictive Error Method (RPEM) 
works extremely well as a stand alone algorithm or as a parameter 
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start-up routine for a block-recursive maximum likelihood (ML) 
algorithm. The block-recursive ML technique gives excellent 
results where several passes can be made over the data block. 

This may not be acceptable for the timing requirements of a real- 
time algorithm. Implementation of a block recursive algorithm 
with a single pass over the data could enhance identification 
speed over pure recursive methods; however, the convergence 
characteristics of such an algorithm require more research. There 
fore, the NRTFA program was configured with RPEM as a stand alone 
algorithm which could be restarted at any point by the operator. 

The equation error, standard ARMA or full ARMAX models can 
be selected by the user (see Appendix A), as described in Section 

2 . 

6.2 THE NRTFA PROGRAM 

The NRFTA subroutine calling structure is given in Figure 

6 - 1 . 

The NRTFA program structure is shown in Figure 6-2. 



Figure 6-1. Near Real-Time Flutter Analysis Calling Sequence. 
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• Form a companion matrix from a(z ) 

• Factor the identified polynomial for 
the discrete roots with the QR 
algorithm 

• Transform discrete roots to the 
continuous domain 

• Sort continuous roots by modulus 

• Compute continuous natural frequency 
and damping ratios 




Figure 6-2. Near Real-Time Flutter Analysis Flowchart 
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Detailed descriptions of the options and input/output argu- 
ments are given with comment statements in the subroutines (see 
Appendix B) ; a brief description of the routines in Fig. 6-1 is: 


NRFTA 


DATAIN 


RPEM 


NSTABL 


RG 


[ HQR3 / 

•< ORHES , 
t ORTRAN 


The main program which allows the user to set 
up the identification problem and calls the 
principal computational and input /output sub- 
routines. 

Opens an input and output data file and calls 
SAVLOD to read the data records. 

Performs an update of the ARMAX model parameters 
and their covariance for one additional input 
and/or output data point. 

Performs a stability check on the C(z poly- 
nomial using the Schur-Cohn test. 

Finds the eigenvalues of a real general matrix, 
used to factor the A(z-l) polynomial by 
operating on the associated companion matrix. 

Supporting eigen-decomposition routines which 
transform a matrix to its Hessenberg form, and 
subsequently to its real triangular form (Schur 
form) with the unsymmetric QR algorithm. 


There are a number of small supporting subroutines which have been 
documented elsewhere [15]. 

There are several default values used to initialize RPEM 
that have been fixed in NRTFA; their values or the way theyare 
initialized could easily be changed in the source code. RPEM 
uses an exponential forgetting factor, discounting past data so 
as to minimize at time t the sum 

j=l 

where A itself is a function of time, 


S£ 



A(j) = (l-a^)A_ + at' A , computed recursively as 
Al A O 

A ( j + 1 ) = (l-a^)A(j) + a^A f . 

The fixed parameters are the initial forgetting factor A (.9), 
the final forgetting factor A^(l), and the rate at which A 

approaches its final value, a. (.97). The values in parenthesis 

A 

are those fixed for flutter application in the early part of 
NRTFA (see the source listings in Appendix B). The above for- 
getting factor parameters achieve a value of .99 (from .9) in 
approximately 75 data points, and .999 in 150 data points. 

In computing the regressors by filtering with the C(z’^) 

polynomial (see Section 3.2.2), a contraction factor is applied 

^ — 1 

to the poles of C(z ), k_c , which reduces the effect of the 
1 ’ f ac ’ 

filtering. k„ is also varied with time by the equation 

I 9.C 


k fac (j+1) = 


%a> c<J) 


+ ( 1 — a. 


L f ac 


where k .. is approaching the value of 1. 

X ctC 

The initial value is set at k„ = .01 with a, = .999. 

f ac k_p 

f ac 

Thus the algorithm starts out as nearly approximately maximum 
likelihood (AML) [7], where the residuals are used directly in 
the regressors, giving good transient convergence characteristics 
and then slowly approaching the RML2 [2] or standard RPEM algo- 
rithm (k„ reaches .9 after approximately 2300 data points), 

X ctC 

for good asymptotic convergence characteristics. 


6.3 MEMORY RESOURCES AND TIMING COMPARISONS 

The active memory requirements for tne NRTFA program and the 
system supporting subroutines on the VAX 11/780 (VMS Version 3.1) 
are shown below in Table 6-1, for the problem definition given in 
Table 6-2. 
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Routine Type 

Memory Requirement (Bytes) 

Main Program 

51,576 

Input /Output Routines 

2,007 

Identification Routines 

2,138 

Eigensystem Routines 

24,138 

System Support Routines 

2,061 

Total Executable Code 

81,920 


Table 6-1. NRTFA Memory Requirements. 

The above memory requirements include storage for the maxi- 
mum sizes shown in Table 6-2 below. 


Parameter 

Maximum Size 

dim. (0) 

60 

Number of modes with an 
A, B, and C polynomial 

10 

Number of modes with 
A and B or A and C 
polynomial 

15 

Number of points per 
data record 

2000 


Table 6-2. NRTFA Program Parameter Limits 
(Effected by Current Array Dimension Bounds) 

To characterize the CPU requirements of different functions 
within the NRTFA program, on the VAX 11/780, the following experi- 
ment was performed. RPEM identifications from 200 data points of 
the A, B, and C polynomials, as well as only the A and B 
polynomial for 2, 4, and 8 modes were performed. The interval 

for outputting the frequency and damping estimates was varied to 
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separate the CPU required for identification from that for factoring 
the polynomial and outputting the desired estimates. In all cases 
the order of the A, B, and C polynomials is twice the number 
of oscillatory modes. The results are tabulated in Table 6-3. 


Number of 
Modes 

Number of 
Points/Batch 

Polynomials 

Identified 

CPU for 
RPEM (sec*) 

CPU for Polynomial 
Factorization and 
, E, output (sec**) 

2 

20 

A , B , C 

1.74/200 

.27/8 

2 

5 

A ,B ,C 

1.74/200 

1.107/32 

4 

20 

A , B , C 

4/200 

2/8 

4 

5 

A , B , C 

4/200 

8/32 

8 

20 

A , B , C 

14.58/200 

11.93/8 

8 

5 

A , B , C 

14.58/200 

47.72/32 

8 

20 

A , B 

7.27/200 

11.93/8 

8 

5 

A , B 

7.27/200 

47.72/32 


Table 6-3. NRTFA Timing Study Results 


Note that the common factor of the time to read the input 
has been subtracted from all cases. These results can be 
much more clearly, graphically. Figure 6-3 shows the CPU 
seconds required for a single data point parameter update 


data 

depicted 
time in 
( 6n 


1.74/200 = 1.74 seconds for 200 updates by the RPEM subroutine. 

.27/8 = .27 seconds for 8 output batch computations. (Factori- 
zation of the polymonial, transformation to continuous fre- 
quencies and damping ratios, as well as output to a terminal 
in the format shown in Appendix A.) 
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parameters of the full ARMAX model for n modes) as a function 
of the number of modes. Similarly the transformation of the 
ARMAX model parameters to continuous frequencies and damping 
ratios is shown in the plot below as a function of the number of 
modes. Clearly overhead becomes more important with a small 
number of modes. 




NO. OF MODES 


Figure 6-3. CPU for Identification and Discrete Polynomial 

Factorization and Output 

(Full ARMAX Model - 6n Parameters for n Modes) 
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These values are consistent with operation counts for these 
two algorithms and the machine cycle time for the VAX 11/780 (VMS 
with a Floating Point Accelerator). The floating point operation 
cycle times for typical machines are shown below in Table 6-4. 


Table 6-4. Millions of Operations Per Second 


Machine 

Vendor Rating 1 (MOPS) 

Performance 

Benchmark^ 

IBM 3033 

4.00 

1.77(D) 3 

Cyber 175 

5.06 

1. 36(S) 

CDC 7600 

10. 

1 . 24 ( S ) 

CDC 6600 

2.5 

. 477 ( S ) 

VAX 11/780 
VMS/FPA 

.831 

. 190 ( S ) 

HP 5451C 
( HP1000 - 
E Series) 

.08 


Apollo 

- 

. 0196 ( S ) 


"Mainframe KOPS Rating," Datamation 1982 

2 

Floating Point (MOPS) for a typical mathematical software 
problem [16] 

3 

(S) or (D) refers to operations in single or double 
precision 

For the 500 sps of the DAST flight data, with, for example, 
a four mode model, and frequency and damping estimates output 
four times a second to the operator, the CPU effort is approxi- 
mately eleven times real time (10 seconds for identification and 
1.2 for factoring the polynomials in one second of flight data). 
These timing results do not include accuracy estimates, which 
require the mode shapes, but beyond this should not require 
significant additional work. The above CPU requirements are 




not as pessimistic as they may first appear for implementing RPEM 
in real time to do flutter analysis. 

Thirty percent of the work in factoring the polynomials can ' 
be eliminated if mode shapes are not desired; however, the 
identification algorithm is where the reduction is needed for a 
real-time implementation. For the modes of interest the sampling 
time could be cut at least in half and with some loss in accuracy 
the sampling rate could be cut by a factor of five. 

The RPEM version implemented here has been a research tool 
that can be streamlined significantly. Even single subscripting 
the arrays could improve efficiency by as much as thirty percent. 
In addition, the implementation of a block recursive algorithm 
could significantly increase the throughput for a given machine 
speed. Such considerations are part of the general analysis 
necessary for actual real-time implementation. 

Other considerations for actual real-time implementation 
include adding a higher level routine to function as an executive 
so that the operator can restart the procedure with function key 
at his console, interfacing a data buffering routine with the 
telemetry system, and displaying the modal parameters and accuracy 
estimates graphically at the operator's console. 

6.4 INSTALLATION INSTRUCTIONS 

The NRTFA program is written in ANSI Standard FORTRAN-77, 
such that conversion to operating systems meeting this standard 
should be straightforward. The storage requirements are modest, 
such that an overlay structure is not required on machines such 
as the Cyber, not having virtual memory. The principal adap- 
tations all have to do with input /output . The subroutine FILES 
which dynamically opens files for the control and output data 
requires a suitably modified OPEN command. The subroutine 
SAVLOD which reads the data files with a specified format into 
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the buffer arrays should be changed to a convenient format for 
the particular installation. The subroutines FILES (not shown as 
it is machine dependent) and SAVLOD are called by the routine 
DATAIN. This file interface is exactly that of MATRIX V [1]. 

The variable EPS in the real general eigensolver (RG) should 
be set to the host machine precision (2^ ^ where t is the 
number of binary bits in the mantissa of a floating point number). 
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SECTION 7 


SUMMARY 


Recent theoretical analysis of recursive identification algo- 
rithms has been applied to flutter test monitoring, and implemented 
in a near-real-time flutter analysis program (NRTFA) shown to be 
highly reliable. This report documents the algorithm and model 
forms used, with the numerical analysis and estimation issues 
considered for its choice. The performance of the recursive 
predictive error method (RPEM) on simulated data compares well 
with predicted convergence rates, accuracy and robustness. The 
results on the DAST flight data show modal estimates that are 
consistent with NASTRAN and ground test vibration results, with 
consideration of the flutter suppression system (FSS) (see 
Gilyard and Edwards [14]). 

Timing studies of the configured FORTRAN program show that 
real-time operation will be possible with careful implementation. 
The development and documentation of this effort has clarified 
the issues to be resolved for a full real-time flutter monitoring 
system. The NRTFA program has been documented sufficiently for 
program modifications, for adaptations to other computers (Section 
6), and for use in current flutter analysis (Appendix A). Sub- 
sequent subsections explain these conclusions more fully. 

7.1 ALGORITHM CHARACTERISTICS 

The RPEM algorithm implemented in the NRTFA program identi- 
fies ARMA-type models (see Section 2.2 for the ARMAX model defi- 
nition), which can be specified for three cases, known inputs, 
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unknown inputs or both. For the first two cases, and sufficiently 
high signal-to-noise ratios in the third case, this algorithm is 
globally convergent. The A(z ^ ) polynomial or equivalently the 
system dynamics is always identified. With known inputs the 
B(z - ^) or equivalently the input distribution matrix is identi- 
fied; with unknown inputs the C(z ) or equivalently the constant 
Kalman gain is identified. 

Two features have been added to effect good transient con- 
vergence behavior. Initially the data is weighted with an 
exponential forgetting factor which discounts past data. This 
factor is changed to approach one, hence tapering off this 
windowing effect after the parameter estimates have stabilized 
somewhat. The NRTFA program also uses a second initial transient 
factor which controls the computation of the estimated gradient. 
This factor constrains the algorithm to be nearly linear regression 
initially [3], [7] which gives a rapid initial convergence rate, 
and then slowly tranforms to a full Gauss-Newton minimization of 
the prediction error [3] for assured asymptotic convergence to 
the global minimum (for the three cases described above). 

7.2 ALGORITHM PERFORMANCE 

Colored sensor noise and higher order dynamics are automat- 
ically accounted for by an overparameterized ARMA-type model. 

The significant conclusion about noise effects is that whenever 
some excitation can be tolerated it is highly desirable for 
improving estimation accuracy. Unknown inputs still yield 
reasonable damping estimates for signal-to-noise ratios above 
approximately ten. Families of curves show the effect of dif- 
ferent turbulence, sensor noise and input levels in Figures 3-2 
and 3-3. The convergence rates and accuracy of estimated param- 
eters from six different simulations compare favorably to pre- 
dicted error bounds in Tables 4-3 through 4-b demonstrating the 
underlying reliability. RPEM is based on a stability analysis of a 
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differential equation associated with the behavior of the param- 
eter estimates as a function of time [6] . Recursive algorithms 
can be made robust to data dropouts and outliers with fairly 
simple logic (see Section 3.3.1), resulting in improved estima- 
tion accuracy (see Table 3-1 for a quantitative comparison). 

For known inputs and a low level of turbulence the recursive 
least squares (RLS - A(z -1 ) and B(z ~^) only) provides very good 
results as long as sufficient modes are modeled to cover unknown 
effects in the vehicle response; for example, with a high-ampli- 
tude swept sine input, a three mode (RLS) model did not give good 
estimates, whereas five and eight mode models did. The unknown 
input results (see Section 5.2) support the general conclusion 
stated previously that a sufficient turbulence level must be pre- 
sent to give a signal-to-noise ratio (SNR) of approximately 10 
for estimation particularity of the damping parameter. With a 
low SNR, even low amplitude known inputs increase the estimation 
of the A(z - ; polynomial and hence frequencies and damping ratio. 
Monitoring of the final portion of the DAST Flight 3 has been 
emulated with the delivered algorithm, showing clearly where the 
flutter mode damping starts its trend to £=-.022 before satura- 
tion of the sensors and controls and eventual failure of the 
right wing. 

7. .3 SOFTWARE IMPLEMENTATION 

Fixed parameters that influence the convergence character- 
istics of RPEM are described in Section 6 such that the NASA 
Dryden personnel can modify the source program (see Appendix B 
for the listings of principal routines) if so desired. Minimal 
necessary installation instructions are also given. Appendix A 
illustrates an example execution of the NRTFA program, with guide- 
lines for input options and description of the modal estimate 
outputs . 
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7.4 RECOMMENDATION FOR FURTHER RESEARCH 


Work is needed in two directions jointly to achieve semi- 
automated monitoring of aircraft flutter in real time. 

1. Further development and testing of algorithms on 
several flight records and simulation runs with a 
variety of errors. Also paramenters like batch 
size, fading factors and model orders need to be 
further investigated. 

2. Development of a real-time flutter analysis hard- 
ware system. 

The hardware system could be used on the ground initially but 
could be extended to be an onboard system as more experience is 
gained with the algorithms and the computation capability advan- 
ces . 

The goal of such a real-time flutter analysis system would 
be to provide a dual environment. The real-time environment mon- 
itors flutter characteristics for safe flight tests. In addition, 
portions of the test are isolated for detailed post-flight analy- 
sis. The second environment is the post-flight processing capa- 
bility. This would be accomplished with a flexible identification 
program to extract more accurate estimates of the frequency and 
damping estimates on critical maneuvers where the real-time 
algorithm indicated very light damping or unusual behavior. Im- 
provements to the real-time algorithms could also be made in 

this "development" or post-flight processing environment, as more 
experience is gained. 

Specific algorithm issues that deserve more attention are a 
better quantification of estimation error effects due to model 
order, more experience with convergence parameter setup values, 
particularly on flight data, the effect and design of prefilters, 
input design for maximum accuracy, and parameter estimation 
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accuracy during the convergence transient phase. 

The timing results of Section 6 indicate that speed needs to 
be a significant concern in the further algorithm analysis, 
quantifying the effects of, for example, not updating the covari- 
ance as often, etc. A comparison of the RPEM, extended Kalman 
Filter (EKF), and maximum likelihood output error formulation 
for a block recursive or semi-batch mode should identify the best 
algorithm for the accuracy, convergence, and speed characteristics 
most suitable for real-time flutter monitoring. 

The joint development of the hardware system with further 
algorithm investigation will be necessary, even though the problem 
appears to be compute bound, because the actual throughput of 
real-time systems can not be fully predicted by study alone. 
Finally, the technology is now available to integrate a fully 
automatic system and validate its performance. 
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APPENDIX A 


USER'S GUIDE 

NEAR REAL-TIME FLUTTER ANALYSIS PROGRAM ( NRTFA ) 


The NRTFA program setup parameters are entered interactively 
as shown below. 

ENTER INPUT FILE NAME... 

U 

ENTER OUTPUT FILE NAME... 

Y 

DO YOU WANT C(Z) POLYNOMIAL, IE. KALMAN GAIN IDENTIFICATION 

(Y OR NO)? 

Y 

DO YOU WANT INPUT POLYNOMIAL B(Z) TO BE IDENTIFIED (Y OR N)? 

Y 

ENTER NUMBER OF DATA POINTS TO BE ANALIZED... 

200 

NUMBER OF MODES TO BE CONSIDERED? 

4 

HOW MANY POINTS PER BATCH? 

20 

The first two prompts are for the excitation and sensor data 
files, where the file names can be up to 32 characters. The input 
file U above contains the sequence of control positions, for 
open loop identification, or the excitation superimposed on a 
closed loop control system for closed loop identification. The 
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excitation channel is always read, whether or not the B(z _1 ) 
polynomial is identified. The output file contains the sequence 
of accelerometer outputs (the present version of RPEM is a SISO 
a I gorithm) . 

The second group of two prompts above defines which two or 
all three of the ARMAX polynomial models are identified (see 
Section 2.2). If there is knowledge of the control channel B(z _1 ) 
should be identified, and if the gust inputs are known to be small 
with wideband accelerometer noise then the recursive least-squares 
(RLS) algorithm (A(z and B(z - '*') identified, without C(z - ^)) 

can be used. In all other cases the full ARMAX model (A(z~^), 

B(z and C(z ^)) should be used. 

The last three prompts of the interactive setup session 
specify the number of data points to be read into buffer arrays, 
the number of oscillatory modes (resulting in 2n or 3n parameters) 
modeled, and how often the frequency and damping estimates are 
output to the terminal. 

The modal estimate output to the terminal resulting from 
I. he session above follows: 
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<• ^ ” ! ' 
r i r i n ; 

i i.* i_ 


S 1 X 7 H H G 1: 

rl 

- c i * 

♦ C 4 * 

5*433 
1 4 ♦ 665 

\J * i w W 

0.356 

* Ur*- * 

* 

7 * a / J 

0*187 

T Lu t 

1 3 , 7 4 9 

0,079 

; ci ; 

0 ♦ \J u 

0,077 

:c 2 i 

9,938 

0,095 

* 



:ci: 

:c4: 

0.493 

13.708 

-0,213 

0.016 

1 C2 5 

♦ 

* 

7.005 

O.OOS 

5C35 

10,038 

0 » 086 

:ci: 

;ct: 

0 ♦ 5od 
13*741 

0,234 

0.071 

* po ♦ 

♦ 

7.151 

0.04 8 

5C35 

10,775 

0.124 

:ci: 

;c4: 

0 . 669 
12.724 

0 ♦ 057 
0*118 

: C2 : 

* 

* 

7.141 

0,095 

5 C 3 5 

9,071 

0.461 

: ci : 
: ca : 

0.617 

13.765 

0,234 

0.092 

: C2 : 
5055 

6,877 

14,519 

0.132 

0 ♦ o wU 

4 0 7* 

* u u * 

* 

• 

9,878 

0,283 

: ci : 

5C45 

0.686 

13.505 

0 . 088 
0.132 

5025 

♦ 

♦ 

6,681 

0,163 

♦ p7 ♦ 
»UU* 

8,927 

0,212 

: ci : 

:C4: 

0.675 

13.763 

0.052 

0,091 

5 025 
5 05 5 

6.570 

14.043 

0.135 

0.218 

♦ p7 * 

♦ L u * 

♦ 

9.557 

0.177 

:ci: 

:C4: 

0.659 

13,414 

0.013 

0,126 

5C25 

♦ 

♦ 

6.483 

0,161 

5C35 

9,111 

0,119 


FORTRAN STOP 

f 

Estimates for the first twenty points are computed before 
output begins, hence these 9 modal output parameter groups. A 
pair of roots are prefixed by whether they are complex (C) or 
two real (R) roots, along with the mode numbers. The modes are 
ordered from smallest to largest by natural frequency. For a 
complex pair the first number is the natural frequency and the 
second the damping ratio. 

For real roots, both are printed. If the discrete roots 
cannot be realized as real continuous modes they are not printed 
at all, for example the second group above, where only two of 
four modes were printed. 
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APPENDIX B 


FORTRAN LISTINGS 


Listings of the major subroutines of the Near Real Time 
Flutter Analysis (NRFTA) Program as described in Section 6 follow, 
ordered according to the calling sequence shown in Figure 6-2. 


PROGRAM RTFA 

DOUBLE PRECISION PMAT ( 60. 60 ) , THETA ( 60 ) . PSID ( 60 ) . FI < 60) . 

1 F'SI ( 60 ) . TSTAB < 21 ) » WORK ( 42 ) » FRF'EM ( 60 > . GRF’EM ( 60 ) . LRF'EM ( 60 ) . 

2 UARRA CIO)? YARRA ( 10) . UCUR. YCUR .LAMBDA . KFACT . LAMDL. ALMDA. 

3 AKFACT .V. RESID.F'INIT .CFACT .PRERR 

c 

c 

C THIS IS THE MAIN PROGRAM FOR REAL-TIME FLUTTER ANALYSIS BY 
C RECURSIVE PREDICTION METHOD (RPEM). 

C HAVING INPUT U(T) AND OUTPUT Y(T) PARAMETER ESTIMATE THETA 
C IS UPDATED EACH TIME SUBROUTINE RPEM IS CALLED. 

C MODEL STRUCTURE USED IS 

C A(Q**-1 )Y(T)=B(Q«-1 )U(T)+C(Q**-1 )E(T) 

C THETA- VECTOR OF ORDER (NARPM+NBRPM+NC) CONTAINING THE 
C . PARAMETER ESTIMATES. 

C THETA=<A(1)»A<2)». . , . A ( NARPM) . B< 1 ) . B<2) . . . . .B(NBRPM) »C(1 ) > 

C C ( 2 ) . . . » > C (NC) ) 

C THIS IS AN INTERACTIVE SEMI-BATCH PROGRAM WITH OPTIONS ON 
C TOTAL NUMBER OF DATA POINTS » NUMBER OF POINTS PER BATCH 

C AND IDENTIFYING [ A < Z ) . B ( Z ) . C (Z ) 3 . CA(Z).B(Z)3 OR CA(Z)f 

C C < Z ) 3 . THE FOLLOWING INPUTS ARE NEEDED. 

C NDATP - NUMBER OF DATA POINTS. 

C NMODE - NUMBER OF MODES TO BE IDENTIFIED. 

C NPRBCH - NUMBER OF POINTS PER BATCH. 

C 

C AFTER NSKIP POINTS. FOR EVERY NPRBCH POINTS NATURAL FREQUENCY 
C AND DAMPING RATIO OF EACH MODE ARE CALCULATED AND PRINTED. 

C THE TRANSFORMED DYNAMIC MATRIX 

C F ( . } =C-A ( ) IN THE FIRST COLUMN . IDENTITY MATRIX! 

C AFTER TRANSFORMING THE EIGENVALUES OF F INTO CONTINOUS 
C DOMAIN THE NATURAL FREQUENCY AND DAMPING RATIO OF COMPLEX 
C ROOTS ARE CALCULATED. IF THE ESTIMATED ROOTS ARE COMPLEX 
C THEN IN OUTPUT ICJ WILL BE PRINTED BEFORE FREQUENCY AND 

C DAMPING RATIO OTHERWISE JRJ AND LOCATION OF REAL ROOTS 
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C WILL BE PRINTED, TO BE ABLE TO PRINT REAL ROOTS THE 

C LINE BELOW LABEL 120 IN THE MAIN PROGRAM 'KDUM=KIMJM3' 

t: SHOULD BE DELETED. 

C IF DAMPING RATIO IS NEGATIVE THAT MODE IS UNSTABLE, 

C 

C F< » ) - DYNAMIC MATRIX 

C FREQ ( )- NATURAL FREQENCY 

DAMP ( )- DAMPING RATIO 
C RLRT( )- REAL ROOTS 

C CLXRT ( )- COMPLEX ROOTS 

C DT - SAMPLING PRIOD 

C 

C INPUT PARAMETERS FOR SUBROUTINES HAVE BEEN DEFINED IN THEM. 

C 

C SUBROUTINE DATAIN READS INPUT U< > AND OUTPUT Y< ), 

C SUBROUTINE RF'EM UPDATES THETA( ) PARAMETER ESTIMATE. 

C SUBROUTINE RG COMPUTES THE EIGEN-VALUES OF F< » ), 

C 

c 

DOUBLE PRECISION F(20» 20) 7 VRG(20»20) >U(500) >Y<500) 

DOUBLE PRECISION FREQ ( 20 > » DAMP (20) » CLXRT <20»2) »RLRT(20) »DT 
DOUBLE PRECISION EPSRG r ER( 20 ) , El < 20 ) > SUPD( 20 ) » DUMP 
INTEGER INPUT < 32 ) > OUTPUT (32) i ID(4) »NR»NC» JOB» IMG 
INTEGER NLOWfNUPf INF0»NV»NARGfTYPE(20)»NDEL(20) 

CHARACTER CHAR ( 10 ) » ANSWR 
C0MPLEX*16 CDUM1 . CDUM2 
DATA IMG»J0B/0?500/ 

DATA NSKIP/40/ 

C 


C EF’SRG= MACHINE EPSILON 

C 

DATA CF ACT/1. OHIO/. IDIM/60/»EPSRG/2.0D-16/ 

DATA NA/20/ »NV/20/ 

DATA DT / , 2292/ 

0PEN<UNIT=5> NAME= ' INPUT ' > STATUS^ 'UNKNOWN ' ) 
0PEN(UNIT=6>NAME='0UTPUT'»STATUS= 'UNKNOWN' ) 

c 

C READING INPUTS 

C 

CALL DATAIN < JOB » IMG > NDATP » U> Y ) 

WRITE ( 6 ? 881 ) 

881 FORMAT ( 2X ? ' DO YOU WANT C(Z) POLYNOMIAL^ IE. KALMAN GAIN ' > 

1 ' IDENTIFICATION (Y OR N)?'> 

READ(5>883)ANSWR 
WRI TE ( 6 » 882 ) 

882 FORMAT ( 2X» 'DO YOU WANT INPUT POLYNOMIAL B(Z) TO BE ' » 

1 ' IDENTIFIED! Y OR N)?') 

READ ( 5 » 883 ) ANSRB 

883 FORMAT (A) 

WRITE ( 61 880) 
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380 FORMAT <2X» ' ENTER NUMBER OF DATA POINTS TO BE ANALIZED . ' ) 
READ(5»*)NDATP 
WRITE ( 6 » 994 ) 

994 FORMAT ( 2X » NUMBER OF MODES TO BE CONSIDERED? . 
READ(5r*)NM0DE 

WRITE ( 6 » 995 ) 

995 FORMAT <2X»' HOW MANY POINTS PER BATCH?') 

READ<5r*.)NPRBCH 

NBACH=NDATP/NF'RBCH 

TW0NMD=2*NM0DE 

C 

c .SET PARAMETERS FOR RF'EM SUB, - 

C 

NARPM-TWONMD 

NBRF'M=TWONMD 

NC=TUONMD 

IF ( ANSWR , NE . ' Y ' ) NC = 0 
IF ( ANSRB »NE * ' Y ' ) NBRPM--0 
ND=0 

MF’=NARPMtNBRPM+NC 

NF’1=NARPM + 1 

NL0W=1 

NUP-NARF’M 

IN I T = i 

IER1=1 

IER2=1 

DO 10 1=1 tMP 
10 THETA ( I ) =0 * ODO 
P I N I T r - 1 • D 4 
LAMBDA-” , 9 DO 
ALMDA= 1 . 0 
LAMDL : *0 , 97 
KFACT =0 ,01 
AKFACT=0*999 , 

I STAB 1 = 1 
ISTAB2=0 
C 

C INITIAL RF'EM CALL 

C 

CALL RF’EM ( THETA ? F’MAT » NARF'M ? NBRF'M » NC > ND ? UCUR ? YCUR > LAMBDA t 

1 KFACT fCFACT> ISTAB1 ? ISTAB2? V » F'RERR rRESID> INIT rF'INIT 

2 iIDIM»IERl»IER2»FI»PSI>PSIDiFRFEMfGRPEM»LRPEM>WQRK> 

3 TSTAB r UARRA , YARRA ) 


WRITE<6»993) 

993 F0RMAT(11X» 'FIRST MODE ' » 1 1 X ? SECOND MODE ' ,1 IX *' THIRD MODE ' i 
1 / 1 1 X » 'FOURTH MODE' » 1 OX » 'FIFTH MODE MIX? SIXTH MODE') 
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c 

c 


MAIN LOOP AND BATCH LOOP SET-UP 


C 

c 

c 


c 


o 


t, 

c 

c 


c 

c 

c 


INIT--0 
I J=0 

DO 30 J=1 jNBACH 
DO 20 1=1 » NPRBCH 
i j=i j+i 
UCUR=U( I J) 

YCUR=Y ( I J) 

IAMBDA=ALMDA*LAMBDA+(1.0D0-ALMDA)*LAMDL 
KFACT=AKFACT#KFACT+( 1 ♦ ODO-AKFACT) 

— RF’EM CALL — 

CALL RPEM( THETA ?F'MAT > NARF'M pNBRF’MpNCpNDjUCURpYCURp LAMBDA » 

1 KFACT f CFACT ? ISTAB1 > ISTAB2? Vf F'RERR»RESID» INIT »F'INIT 

2 ,IDI«»IERl»IER2tFI»PSI»P8ID»FRPEMf0RPEN»LRPEMfH0RK» 

3 TSTABf UARRAf YARRA) 

20 CONTINUE 


TRANSFORMED DYNAMIC MATRIX SET-UP — 

IF( I J*LT . NSKIP) GO TO 30 
DO 50 IDUM=1 pNARF'H 
DO 50 JDUM=1 jNARF'M 
F ( IDUM • JDUM ) =0 * ODO 
IF(. JDUMtEO* IDUM+1 ) F < IDUM > JDUM ) =1 . 0 
50 CONTINUE 

DO 40 K=1 pNARF'M 
F (K»1)=-THETA(K) 

40 CONTINUE 

CALL RG TO COMPUTE EIGENVALUES 

CALL RG( F ?NA» NARPMp VRG» NVjNLOUpNUP » ER»EI » TYPEpSUF’Dp 
1 NDELpNF'IpNBLOCKpEF'SRGpINFO) 

TRANSFORMING ROOTS TO CONTINOUS DOMAIN - 


K = 0 

KDUM=0 
KDUM1 = 0 
60 K=K+ 1 

IF(K.GE.NARPM) GO TO 90 
KP1=K+1 

IF(ABS(F<K»KPl)),LT.1.0D-07) GO TO 110 
KDUM=KDUM+1 

CDUM1=DCMF’LX<F{KpK) * F < K f KP 1 ) ) 
CDUM2=CDL0G ( CDUM1 ) /DT 
CLXRT(KDUM,1)=DREAL(CDUM2) 

CLXRT (KDUM»2)=DIMAG(CDUM2) 

K=K+1 
60 TO 60 
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110 I F < ABS ( F ( K » K ) ) ♦ GE * 1 ♦ ODO ) GO TO 60 
IF(F (K»K) . LT.O.OBO) GO TO 150 
KDUM1=KDUH1+1 

RLRT (KDUH1 )=BLOG(F (K»K) )/BT 
GO TO 60 
150 KDUM=KDUH+1 

CLXRT ( K'BUM » 1 ) =DLGG ( -F ( K » K > )/DT 


CLXRT (KBUM»2)=3. 1415/DT 
GO TO 60 
90 CONTINUE 

BO 80 K=1 »KDUM 

IF ( ABS ( CLXRT (K» 2) ) «GT . 1 .0D-07 ) GO TO 70 
FREQ'(K) =CLXRT (K » 1 ) 

DAMPC K)=CLXRT (K» 1 ) 

CHAR ( K ) = 7 R ' 

GO TO 80 

70 FREQ(K)=DSQRT (CLXRT <K, 1 >*CLXRT <K> 1 ) +CLXRT1 K,2> #CLXRT <K»2> ) 
DAMF’(K)=-CLXRT (Kr 1 )/FREQ(K) 

CHAR <K) = 'C ' 

80 CONTINUE 

SORTING THE ROOTS 

KDUM3=KDUf1 
BO 140 K=1 >KDUM-1 
BO 140 K F‘ 1 = K + 1 »KBUM 

IF ( ABS (FREQ(K) -FREQ (KP1 ) ) ♦ LT «0«0001 )KDUM3=KBUM3-1 
IF ( ABS (FREQ( K) -FREQ (KP1) ) ♦ LT ♦ 0.0001 )FREQ(KP1 ) =1 . 0D+10 
IF (FREQ(K) *LT *FREQ<KP1 ) ) GO TO 140 
BUMP=FREQ(K> 

FREQ ( K ) =FREQ ( KP1 > 

FREQ(KP1)=DUMP 
BUMF'=DAMP ( K ) 

BAMPCK)=BAMP(KP1> 

DAMP ( KF’l ) =DUMP 
ANSWR=CHAR ( K ) 

CHAR ( K ) =CHAR ( KP1 ) 

CHAR1KP1 ) =ANSWR 
140 CONTINUE 
KDUH=KBUM3 

IF(KBUH.GT.NMOBE) GO TO 120 
DO 130 K=1 i KDUM1 
KDUM-KDUM+1 
FREQ (K BUM )=RLRT (K) 

DAMP(KBUH)=RLRT(K+1 ) 

I F < K . GT ♦ KDUM1 ) DAMP (KDUM ) =0 , 0 
130 CHAR(KDUM)= , R' 

120 IF (NDUM.GT *NMOBE) KBUM=NMODE 
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OUTPUT PRINT-OUT 


C 

c 

KDUM=KDUH3 
WRITE ( 6.883) 

WRITE<6f992)!CHAR<K>»K»FREQ<K)»DAHP<K')»K*l»KDUM> 
992 FORMAT ! 3(3X. ' ♦ ' » A1 » 1 1 » ' * ' . 2F8.3 )) 

30 CONTINUE 
STOP 
END 


SUBROUTINE DATAIN ( JOB > IMG » NBATP » U t Y ) 

C 

C SUBROTINE DATAIN 

C 

DOUBLE PRECISION Y<1>.U!1) 

INTEGER ID ( 4 ) . NR » NC r JOB. IMG 
INTEGER INPUT (32) > OUTPUT (32) 

C 

C INPUT FILE NAMES 

C 

WRITE (6» 9000) 

9000 FORMAT (/»' ENTER INPUT FILE NAME . . . ' ) 
READ(5»9010»END=2000) INPUT 
, 9010 FORMAT (32A1) 

WRITE < 6^9030 ) 

9030 FORMAT (/» ' ENTER OUTPUT FILE NAME,..') 
READ (5. 9010. END = 20 00 ) OUTPUT 

C 

C ASSIGN FILES 

C 

CALL FILES! 2. INPUT ) 

CALL FILES! 1 1 .OUTPUT ) 

C 

C — LOAD MATRIXX FILE 

C 

CALL SAULOEK 2 » ID. 500 r NR. NC? IMAG . JOB » U . U) 
CALL SAVL00<11.IB.500>NR.NC.IMAG.J0B.Y.Y> 
NDATP=NR 

C 

2000 CONTINUE 
C 

C CLOSE FILES 


CALL FILES! -2 . INPUT ) 

CALL FILES! -1 1 r OUTPUT) 

RETURN 

END 

SUBROUTINE SAOLOD ( LUNIT . ID. HA . M , N . IMG. JOB* XREAL . XIMAG ) 
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c 

c 


SUBROUTINE SAVLQAD 


INTEGER LUNIT>ID<4)»MiN»IMG> JOB 
DOUBLE PRECISION XREALI 1 ) > XIMAG( 1 ) 

C 

C IMPLEMENT SAVE AND LOAD 

C LUNIT = LOGICAL UNIT NUMBER 

C ID = NAME > FORMAT 4A1 

C M, N = DIMENSIONS 

C IMG = NONZERO IF XI MAG IS NONZERO 

C JOB = 0 FOR SAVE 

C = SPACE AVAILABLE FOR LOAD 

C XREAL i XIMAG = REAL AND OPTIONAL IMAGINARY PARTS 

C 

C SYSTEM DEPENDENT FORMATS 

101 FORMAT ( 4A1 » 3 1 4 ) 

102 FORMAT ( 4Z18 ) 

IF (JOB ,6T, 0) GO TO 20 

c; 

C SAVE 

C 

10 WRITE ( LUNIT » 101 ) ID»M»N»IMG 
DO 15 J = If N 
K = < J-l >*MA + 1 


L - K + M - 1 

WRITE ( LUNIT » 102) ( XREAL ( I ) » I=K > L ) 

IF (IMG .NE. 0) WRITE(LUNIT » 102) (XIMAG(I) »I=KfL) 

15 CONTINUE 
RETURN 

C 

C LOAD 

C 

20 READ (LUNIT » 101 jEND=30»ERR=29) ID,H>N>IMG 
IF (M*N . GT . JOB) GO TO 30 
DO 25 J = 1, N 
K = <J-1)*M+1 
L = J*M 

READ(LUNIT > 102»END=30) (XREAL ( I ) , I=K» L ) 

IF (IMG .NE. 0) READ (LUNIT i 102»END=30) ( XIMAG ( I ) , I=K > L ) 
25 CONTINUE 
RETURN 

29 WRITE(6r*) 'ERROR IN READING FILE' 

C . 

C END OF FILE 

C 

30 M = 0 
N = 0 
RETURN 
END 


115 



SUBROUTINE RPEM(THETA*F’*NAiNB*NC*ND*U* Y*LAMBDA*K*C> 

C ISTAB1 * ISTAB2* V»EF’S» EF'Sl » INIT* PQ* IDIM* IERI j IER2* 

C FI jF'SI >PSID*F*G*L*WORK»TSTAB* UARRA* YARRA) 

RPEM SUBROUTINE 

RECURSIVE PREDICTION ERROR METHOD 

THIS IS A MODIFIED VERSION OF THE ROUTINE ORIGINALLY 
DEVELOPED BY T.SODERSTROM. THE MODIFIED FORM HANDLES 
VARIABLE SIZES OF POLYNOMIALS A * B » AND C * THE ROUTINE 
PROVIDES VARIOUS OPTIONS AS DESCRIBED BELOW, 

THE SUBROUTINE PERFORMS THE MODIFICATION OF THE 
PARAMETER ESTIMATES THETA FOR ONE SAMPLING INTERVAL 
A NEW CALL TO RPEM MUST BE HADE FOR EVERY NEW 
SAMPLING INTERVAL 
MODEL STRUCTURE USED IS 

A(GM-1 )Y(T)=B(Q**-1 )U(T-NB)+C(Q*#-1 )E(T> 

THETA - VECTOR OF ORDER (NA+NB+NC) CONTAINING THE 
PARAMETER ESTIMATES, 

THETA=(A(i ) » . . . »A(NA) *B(1) * . . , *B(NB) *C(1) * , , . . *C(NC) > 
THETA IS CHANGED IN THE SUBROUTINE 
P - SYMMETRIC MATRIX OF ORDER! NA+NB+NC ) 

P IS USED IN THE U-D FORM, (THIS U IS DIFFERENT 
FROM THE INPUT VARIABLE U(T).THE ARGUMENT OF THE 
ROUTINE CONTAINS U<T>>. 

P=U#D#U (TRANSPOSED) WITH D DIAGONAL AND U UPPERY 
TRIANGULAR, THE ELEMENTS OF D ARE STORED IN THE 
DIAGONAL OF P. THE ELEMENTS OF U ARE STORED IN THE 
UPPER TRIANGULAR PART OF P, P IS CHANGED IN THE 
SUBROUTINE, 

U -THE LAST INPUT VALUE 

Y - THE LAST OUTPUT VALUE 

LAMBDA - THE FORGETTING FACTOR (TO BE ENTERED) 

O.LT, LAMBDA, LE.l , 

K - THE CONTRACTION FACTOR 
0 »LT , K*LE * 1 , 

THIS FACTOR CONTRACTS THE ROOTS OF THE 
C - POLYNOMIAL , (THIS IS TO BE ENTERED) 

THIS FACTOR IS USED IN FILTERING OF THE DATA 
C - PARAMETER USED FOR THE REGULARIZATION (SHOULD BE 
CHOSEN RATHER LARGE, TO BE ENTERED) 

THIS LIMITS THE MAXIMUM VALUE OF THE DIAGONAL 
ELEMENTS OF THE P - MATRIX 

ISTAB1 - FLAG (TO BE ENTERED) FOR STABILITY TEST OF C(Z) 
IF ISTAB1=0* NO MONITORING (STABILITY TEST AND 
STEP SIZE REDUCTION) IS PERFORMED, 

IF ISTAB1 ,NE.O» MONITORING IS PERFORMED 
ISTAB2 - INTEGER AT RETURN GIVING THE NUMBER OF STEP 
SIZE REDUCTIONS PERFORMED 

V - LOSS FUNCTION- SUM OF SQUARED PREDICTION ERRORS, 

TRANSIENT PHASE IS INCLUDED 



c 

c 

I'. 

c 

c 

c 

c 

c 

c 

c 

c 

c 


EPS - GIVES THE PREDICTION ERROR ON RETURN 
EF’Sl - GIVES THE RESIDUAL ON RETURN PROVIDED UPDATE USES 
RESIDUALS BY IER1=1 

I N IT - FLAG TO BE USED FOR STARTING THE RECURSION, IF 

INIT=0» ALL PARAMETERS ARE UPDATED. IF INIT.NE.O, 
APPROPRIATE INITIAL VALUES ARE FIRST SET AND THE 
PARAMETERS ARE UPDATED USING THE AVAILABLE DATA 
U,Y 

PO - SCALAR PARAMETER USED TO GIVE AN INITIAL VALUEUO 
BE ENTERED WHEN INIT.NE.O) 

IF INIT.NE.O, P=PO*I 
IDIM - DIMENSION PARAMETER 


C IER1 FLAG TO BE ENTERED, IF IER1=0, PREDICTION ERRORS 

C ARE USED IN PLACE OF THE RESIDUALS. IF IER1.NE.O, 

C THE ALGORITHM USES THE RESIDUALS. 

C IER2 - FLAG TO BE ENTERED. IF IER2=0, FILTERING IS 

C NOT USED IN THE ALGORITHM. IF IER2.NE.0, 

C FILTERING WITH POLYNOMIAL C(Z> IS PERFORMED, 

0 ABS(ND) ,LT . 10 

C IF IER1 =0 AND IER2=0, THE ALGORITHM REDUCES TO THE 

C EXTENDED LEAST SQUARES CASE. 


C IF IER1=0» AND NB=NC=0, THE SUBROUTINE REDUCES TO A 

C SIMPLE RECURSIVE LEAST SQUARES ESTIMATE OF AR MODEL, 

C 

C IF IER2=0» THE ALGORITHM PERFORMS AS RML1 

DOUBLE PRECISION THETA ( 1 ) »P(IDIM , 1 )» PSID (1 ) , 

C FI(l>,PSI(i),TSTAB<l),WORKU),F(l>,G(l>»Ui,Yl»El,AL 
C * Cl , EPS , EPSlf ALFA, PO, V, C, U,Y,AMY,DD, S,BETA, GAMMA 
DOUBLE PRECISION LAMBDA, K , L< 1 ) * UARRAI 1 ) , YARRA< 1 > 
NN=NA+NB+NC 

C TEST FOR INITIALIZATION 

IF(INIT.EQ.O) GO TO 100 


0=0 , DO 

DO 10 1=1, NN 
DO 10 J=i,NN 
10 P( I , J) =0 » DO 

DO 15 1=1,10 
Y ARRA ( I ) =0 . DO 
15 UARRAI I >=0.B0 

DO 20 1=1, NN 
P < I , I ) =P0 
L ( I ) =0 » DO 
f I ( I ) =0 , DO 


20 

PSI ( I )=0.D0 


NS=3*NC 


DO 30 1=1, NS 

30 

PSIDI I )=0.D0 


RETURN 

C 
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100 NDS-ND 

IF < NOS ,LT. 0) 60 TO 103 
IF (NH *LT . 1) GO TO 102 
DO 101 1 = 1 » ND 
J=NDS+2-I 

101 UARRA( J)=UARRA( J-l ) 

102 UARRA ( 1 ) =U 
U=UARRA(ND+1 ) 

GO TO 108 

103 ND5=-NB8 

DO 105 I=1tNDS 
J=NDS+2-I 

105 YARRA( J)=YARRA< J-l ) 

YARRA ( 1 ) = Y 
Y=YARRA(NDS+1 ) 

108 CONTINUE 

COMPUTE PREDICTION ERROR 
EPS=Y 

DO 110 1 = 1 t NN 

110 EPS=EPS-FI(I)ITHETA(I) 

C COMPUTE NEW PARAMETER ESTIMATES 

AMY = l t 

C TEST FOR NEED OF MONITORING 

IF < NC . EG * 0 ) GO TO 200 
IF ( ISTAB1 *EQ .0 ) GO TO 200 


ISTAB2=0 

120 DO 130 1=1 >NC 
NI=NA+NB+I 

130 T8TAB<m>=THETA(NI)+L(NI>*EPS#AMY 

TSTAB ( 1 ) = 1 * 

C TEST FOR STABILITY OF C(Z) 

CALL N3TABL < TSTAB j NC » WORK * 1ST ) 

IF ( 1ST .EQ.0)G0 TO 200 
AMY = AMY/2 * 

ISTAB2=ISTAB2+1 

IF ( I3TAB2 .GT . 25) RETURN 

GO TO 120 

C UPDATE PARAMETER ESTIMATES 

200 DO 210 1 = 1 f NN 

210 THETA < I ) =THETA ( I ) +L< I ) *EPS*AMY 

I F ( IER1 . EG * 0)G0 TO 250 
C COMPUTE RESIDUALS 

EPS1=Y 

DO 220 1 = 1 ? N M 

220 EPS1=EPS1-FI(I)*THETA(I) 

250 IF ( IER1 .EG. 0 )EF'S1=EPS 
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c COMPUTE FILTERED SIGNALS AND STORE IN PSID ARRAY 

Y1=Y 

U1=U 

E1=EPS1 

IF ( IER2.EG.0)G0 TO 670 
IF(NC.EG.O)GO TO 670 
DO 620 1=1 »NC 
CI=THETA(NA+NB+I)*K**I 
Y1»Y1+CI*PSID< I > 

U1=U1-CI#PSID(NC+I) 

620 E1=E1-CI*PSID(2*NC+I> 

IF < NC . EG . 1 ) GO TO 650 
C UPDATE PSID VECTOR 

DO 630 I - 2 r N C 
1 1 =NC+2-I 

PSID < 1 1 ) =PSID( 1 1-1 > 

I2=NC+NC+2-I 
PSID(I2)= PSID (12-1) 

I 3=NC+NC+NC+2-I 
630 P8ID( I3)=PSID( 13-1 > 

650 PSID(1)=-Y1 

PSID(NC+1 )=U1 
F’SID(NC+NC+1 )=E1 

C UPDATE VECTORS FI AND PSI 

670 CONTINUE 

IF ( NA .EG . 1 ) GO TO 720 


DO 700 J=2 1 NA 
I =NA+2~ J 
FI ( I > =FI ( 1-1 ) 

700 PSI ( I ) =PSI ( 1-1 ) 

720 FI ( 1 ) =-Y 

F'SI(i)=-Yi 

IF ( NB .EQ .0)60 TO 750 
IF ( NB .EG » 1 )G0 TO 740 
DO 730 J=2 1 NB 
I =NA+NB+2- J 
FI < I )=FI ( 1-1 ) 

730 PSI ( I ) =PSI ( 1-1 ) 

740 FI ( NA+1 ) =U 

PSI ( NA + 1 ) =U1 


750 IF (NC. EQ . 0) GO TO 780 

IF ( NC » EG • 1 ) GO TO 770 
DO 760 J=2 ? NC 
I =NA+NB+NC+2- J 
F I < I ) =FI ( 1-1 ) 

760 PSI ( I ) =F'SI < I — 1 ) 

770 FI ( NA+NB+1 ) =EPS1 

PSI ( NA+NB+1 )=E1 



780 

C 

800 

810 

C 


820 

830 

C. 

835 

840 

t 


CONTINUE 

COMPUTE GAIN VECTOR It UPDATE P AND V 
DO 810 I=2»NN 
J=NN+2-I 
ALFA=PSI ( J ) 

J1=J-1 

DO 800 KK = 1 » J 1 
ALFA=ALFA+P(KK»J)*PSI(KK> 

F < J)=ALFA 
G< J)=P( Jr J)*ALFA 
G < 1 ) — P < 1 » 1 )#PSI < 1 ) 

F(1)=PSIU> 

ALFA=L.AMBDA+F { 1 )*G( 1 ) 

GAMMA=0 * 

IF ( ALFA. GT,0,)GAMMA=1. /ALFA 
I F ( G ( 1 ) « NE ♦ 0 » ) P ( 1 » 1 ) =GAMMA#P < 1 » 1 ) 

DO 830 J=2fNN 
BETA=ALFA 
DD=G( J) 

ALFA=ALFA+DD*F<J> 

IF < ALFA *EQ *0 * ) GO TO 835 
AL=-F< J))KGAMMA 
J 1 = J - 1 

DO 82<' $ = 1 > J1 
S=P(I» J) 

P<I> J)=S+AL*G(I> 

G(I>=G(I)+DD*S 
GAMMA=1 ./ALFA 

P ( J ? J ) =BET AfcGAMMAfcF’ ( J » J ) /LAMBDA 
P(J,J)=MIN(P(JiJ>»C) 

CONTINUE 

V=V+EPS**2/ALFA 

CONTINUE 

DO 840 1 = 1 > NN 

L ( I )=G( I ) /ALFA 

RETURN 

END 


120 



SUBROUTINE NSTABL ( A> N » U » 1ST ) 

TEST FOR STABILITY 

DOUBLE PRECISION A(1)»W<1)»AL 


IST = 1 
N1=N+1 
DO 1 I=1»N1 
W( I )=A< I) 

W ( Nl + I ) =0 1 
K>0 

IF(K«EQ.N)G0 TO 99 

NK1=N-K+1 

DO 11 J=1»NK1 

W( N1+ J) =W (NK1-J+1 ) 

IF(W(N1tNK1 ) .EQ»0* )G0 TO 98 

AL=U(NK1)/W(N1+NK1) 

IF( ABS( AL) *GE . 1 «0)G0 TO 98 
MK=N-K 

DO 12 J = 1 » NK 
W< J)=W(J)-AL*W(N1+J> 

K=K+1 
GO TO 10 
RETURN 
I ST = 0 
RETURN 
END 



SUBROUTINE RG ( A * NA > N * V > NV * NLOW * NUF' * ER * El > T YF'E * SUF'D * 

RG SUBROUTINE 

1 NHEL » N P 1 » NBLOCK * EPS* INFO ) 

INTEGER NA* N*NV* NLOW* NUP* NBLOCK* TYPE (N) *NDEL(NP1) 

DOUBLE PRECISION EPS* A(NA* N> * V(NV*N) >ER(N) *EI <N) »SUPD(N) 

RG FINDS THE EIGEN VALUES OF A MATRIX. ALSO IT HAS CAPABILITY 
OF FINDING EIGEN VECTORS TOO. TO BE ABLE TO FIND EIGEN- 
VECTORS SUBROUTINE ELMBTR WHICH HAS BEEN COMMENTED OUT MUST 
BE CALLED. 

ON ENTRY 

A CONTAINS THE INPUT MATRIX 

NA*N LEADING DIMENSION AND COLUMN DIMENSION OF A 

NLOW* LOW AND HIGH INDEX OF DIAGONAL POSITIONS OF THE 

NUP PORTION OF A TO BE REDUCED TO SCHUR FORM 


ON RETURN 

A DESTROYED ON OUTPUT 

V CONTAINS THE PRINCIPAL VECTORS OF A 

NV LEADING DIMENSION OF V 

ER * El REAL AND IMAGINARY PARTS OF THE N EIGENVALUES OF A 

TYPE (INTEGER) FLAG ARRAY INDICATING THE TYPE OF EIGENVALUE 
(0 => REAL * *1 => POS. CONJUG* 2 = > NEG CONJUG) 

SUPD SUPERDIAGONAL COUPLING ELEMENTS (INDEXED ON THE COLUMN 
POSITION) 

NDEL FLAG ARRAY RECORDING THE LOCATION AND SIZE OF THE 

BLK TRIANGLES* (NDELd + l) IS THE BOTTOM CORNER OF BLK I 


NBLOCK NUMBER OF JORDAN MATRICES FOUND 



C INFO = 0 NORMAL RETURN 

C 

C =. 1 LACK OF CONVERGENCE IN INITIAL GR DECOMPOSITION 

C TYPE ( I ) = -1 ARE THE UNCONVERGED EIGENVALUES 

C 

C EPS IS THE RELATIVE MACHINE PRECISION* I.E. 2**<1-T> 

C WHERE T IS NUMBER OF BINARY DIGITS OF FLOATING F’T 

C MANTISSA < FOR IBM 16**(l~i4> = 2.22B-16) 

C 


I 0RDER=0 
FNORM=O.ODO 
DO 10 1 = 1 > NA 
DO 10 J — 1 » N 

FNORM=FNORM+A (I * J ) *A ( I * J ) 


10 CONTINUE 

CALL ORTHES ( NA ? N , NLOW , NUF' , A , SUPD ) 

CALL QRTR AN ( NA , NV , N , NL OW , NUF' » A * SUPD , V > 

CALL HQR3( A,NA,N,NLQW,NUF', IORDER, V,NV»ER»EI , TYPE, SUF’D, EPS) 

C 

C , . ♦ SET UP FLAGS FOR ELMBTR 
C 

NDEL ( 1 ) =0 

NBL0CK=1 

DO 100 1 = 1, N 

SUPD ( I ) =0 ♦ ODO 

IF ( TYPE ( I ) »LT . 0 ) INFO=I 

IF ( TYPE ( I ) » EG . 2 ) GO TO 100 

NBL0CK=NBLQCK+1 

NDEL ( NBLOCK ) =NDEL( NBLOCK- 1 ) +TYPE ( I )+l 
100 CONTINUE 

NBL0CK=NBL0CK-1 

C 

C , ,, ELIMINATE TO BLOCK TRIANGLAR FORM 
C . * TO SAVE SOME CPU TIME.,,. 

C 

C CALL ELMBTR<A,NA»N,V,NV,ER,EI,TYPE»SUPD»NDEL,NP1» 

C 1 NBLOCK, FNORM, EPS) 

RETURN 

END 
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